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Simulation  of  Integrated  Surface-Water/Ground- 
Water  Flow  and  Salinity  for  a  Coastal  Wetland  and 
Adjacent  Estuary 

By  Christian  D.  Langevin,  Eric  D.  Swain,  and  Melinda  A.  Wolfert 


Abstract 


The  SWIFT2D  surface-water  flow  and  transport  code,  which  solves  the  St.  Venant  equations  in  two  dimensions, 
was  coupled  with  the  SEAWAT  variable-density  ground- water  code  to  represent  hydrologic  processes  in  coastal  wet¬ 
lands  and  adjacent  estuaries.  The  integrated  code  was  applied  to  the  southern  Everglades  of  Florida  to  quantify  flow 
and  salinity  patterns  and  to  evaluate  effects  of  hydrologic  processes.  Results  indicate  that  most  surface  water  within 
Taylor  Slough  flows  through  Joe  Bay  and  into  Florida  Bay  through  Trout  Creek.  Overtopping  of  the  Buttonwood 
Embankment,  a  narrow  but  continuous  ridge  that  separates  the  coastal  wetlands  from  Florida  Bay,  does  occur  in 
response  to  tropical  storms,  but  the  net  overflow  is  only  1.5  percent  of  creek  discharge.  The  net  leakage  rate  for  the 
coastal  wetland  is  about  zero  with  nearly  equal  upward  (17.1  cm/yr)  and  downward  (17.4  cm/yr)  rates.  During  the  dry 
season,  the  coastal  wetland  increases  in  salinity  to  30-35  practical  salinity  units  but  is  flushed  each  year  with  the  onset 
of  the  wet  season.  Model  results  demonstrate  that  surface-water/ground-water  interactions,  density-dependent  flow, 
and  wind  affect  flow  and  salinity  patterns. 


1.  Introduction 

Coastal  wetlands  are  a  difficult  hydrologic  environment  to  represent  with  a  numerical  model  because  of  the  large 
number  of  contributing  hydrologic  processes,  shallow  hydraulic  gradients,  and  variable-density  flow  conditions. 
Existing  numerical  modeling  strategies  have  been  developed  for  either  the  freshwater  wetland  system  or  the  estuary, 
but  simulations  rarely  span  both  domains.  Recently,  distributed-parameter  physics-based  computer  codes  have  been 
developed  to  simulate  coupled  surface-water  and  ground- water  flow  for  inland  freshwater  systems.  Examples  include: 
InHM  (VanderKwaak,  1999;  VanderKwaak  and  Loague,  2001),  MIKE  SHE  (Graham  and  Refsgaard,  2001), 
MODHMS  (HydroGeoLogic  Inc.,  2003;  Panday  and  Huyakorn,  2004),  and  WASH123  (Yeh  and  Huang,  2003).  To 
simplify  surface  and  subsurface  coupling  techniques  and  to  minimize  computer  runtimes,  many  integrated  models  use 
the  diffusive  wave  approximation  to  the  St.  Venant  equation  to  represent  overland  flow.  The  diffusive  wave  approx¬ 
imation,  in  which  the  convective  and  local  acceleration  terms  are  neglected,  is  normally  a  valid  approximation  for 
inland  systems  due  to  relatively  high  frictional  resistances,  small  flow  velocities,  and  shallow  flow  depths.  Most  inte¬ 
grated  models  are  also  based  on  the  assumption  of  constant  fluid  density,  and  thus  their  applicability  to  coastal  regions 
is  questionable  unless  it  can  somehow  be  shown  that  model  results  are  insensitive  to  density  variations.  Conversely, 
estuary  and  oceanic  models  typically  solve  the  full  St.  Venant  equations  because  the  convective  and  local  acceleration 
terms  are  significant  under  tidal  and  wind-driven  conditions.  Furthermore,  most  estuarine  and  oceanic  models  contain 
options  for  including  the  effects  of  density  on  surface-water  flow,  and  have  transport  capabilities  to  simulate  salinity. 
Estuarine  and  oceanic  models,  however,  normally  assume  ground-water  exchanges  are  negligible,  or  that  the 
exchanges  can  be  represented  as  a  simple  source  term  (Wang  and  others,  2003;  Brown  and  others,  2003).  Thus,  most 
of  the  existing  codes  are  not  well  suited  to  represent  both  the  inland  and  marine  systems,  and  the  coastal  wetlands  that 
separate  them. 
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Land 

surface 


This  paper  describes  the  development  and  application 
of  an  integrated  surface-water/ground-water  flow  and  sol¬ 
ute-transport  code  designed  to  simulate  two-dimensional 
overland  flow  and  three-dimensional  fully  saturated 
ground- water  flow.  The  integrated  code  was  designed 
specifically  for  the  coastal  wetland  transition  zone 
between  inland  freshwater  systems  and  marine  systems 
(fig.  1).  Surface-water  flow  and  transport  are  simulated 
using  the  SWIFT2D  two-dimensional,  finite-difference 
hydrodynamic  code  originally  designed  for  estuaries 
(Leendertse,  1987).  The  SEAWAT  three-dimensional, 
finite-difference  code  is  used  to  simulate  variable-density 
ground- water  flow  (Guo  and  Langevin,  2002).  The  two 
models  are  explicitly  coupled  with  a  one-timestep  lag 
using  a  variable-density  form  of  Darcy’s  law  for  flow 
exchange  and  non-diffusive  salt  flux  between  models. 
The  paper  first  describes  the  governing  equations  for  flow 
and  transport  in  both  systems  and  then  presents  the  numer¬ 
ical  procedure  for  implementing  the  two  codes  in  a  coupled  framework.  Lastly,  the  integrated  code  is  applied  to  the  southern  Ever¬ 
glades  of  Florida  and  northeastern  Florida  Bay  to  quantify  flow  and  salinity  patterns  for  a  7-year  period  (1996-2002)  and  to  examine 
the  effects  of  selected  hydrologic  processes. 
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Figure  1.  Conceptual  model  for  surface-water  and  ground-water 
interactions. 


2.  Governing  Equations 

The  subsequent  governing  equations  are  well  described  in  the  literature,  and  have  been  selected  to  represent  hydrologic  pro¬ 
cesses  in  coastal  wetlands  and  adjacent  estuaries.  The  two-dimensional  vertically  averaged  flow  equations  are  used  for  the  surface 
flow  as  a  compromise  that  allows  better  horizontal  resolution  at  the  cost  of  vertical  resolution.  This  is  justified  by  the  observation 
that  in  coastal  wetlands,  it  is  important  to  accurately  represent  topographic  relief,  because  variations  in  ground- surface  elevations 
are  of  the  same  order  as  water  depths,  while  the  shallow  depths  make  baroclinic  driving — a  main  cause  of  third-dimension 
flow — highly  ineffectual.  The  equations  used  to  couple  the  surface-water  model  with  the  ground- water  model  assume  that  unsat¬ 
urated  zones  are  thin  to  absent,  and  leakage  to  the  water  table  can,  therefore,  be  treated  as  instantaneous.  This  assumption  may  limit 
the  approach  to  areas  with  shallow  water  tables  and  highly  porous  materials. 


2.1  Surface-Water  Flow  and  Solute  Transport 

The  governing  equations  for  a  shallow  surface-water  system  consist  of  conservation  of  mass,  volume,  and  momentum.  Leen¬ 
dertse  and  Gritton  (1971)  and  Leendertse  (1987)  present  the  following  governing  equations,  which  were  modified  by  Swain  and 
others  (2004)  to  include  aerially  distributed  sources  and  sinks,  describing  the:  (1)  conservation  of  water  volume,  (2)  conservation 
of  momentum  in  the  v-direction,  (3)  conservation  of  momentum  in  the  y-direction,  and  (4)  solute  mass  transport: 
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where  h  is  water  stage  [L],  d  is  water  depth  [L],  vx  and  vy  are  vertically  averaged  velocities  in  the  v-  and  y-directions  [LT1],  qsg  is 
a  source/sink  term  representing  the  volumetric  exchange  between  surface  water  and  ground  water  per  unit  area  [LT1],  qr  is  a  rainfall 
source  term  representing  the  volumetric  rate  per  unit  area  [LT1],  qet  is  an  evapotranspiration  sink  term  representing  the  volumetric 
rate  per  unit  area  [LTl],fis  the  Coriolis  parameter  [T1],  g  is  gravitational  acceleration  [LT2],  p  is  water  density  [ML-3],  R  is  the 
bottom-stress  coefficient  [T1],  Cj is  the  wind-stress  coefficient  [L°],  pa  is  air  density  [ML-3],  W is  wind  speed  [LT1],  \|/is  the  angle 
between  wind  direction  and  the  positive  y-axis  [degrees],  k  is  the  horizontal  momentum-exchange  coefficient  [L2T^,  C  is  solute 
concentration  for  a  conservative  non-reactive  constituent  [ML-3],  Dx  and  Dy  are  the  dispersion  coefficients  in  the  v-  and  y-directions 
[L2T!],  Csg  is  the  leakage  concentration  between  surface  water  and  ground  water  [ML-3],  and  Cr  is  the  solute  concentration  of  rain¬ 
fall.  In  this  paper,  the  source  concentration  for  rainfall  and  the  sink  concentration  for  evapotranspiration  are  both  assumed  to  be 
zero,  because  C  represents  salinity  concentration,  which  is  considered  conservative  and  non-reacting.  The  transport  equation  (eq.  4) 
can  easily  be  extended  to  represent  reactive  and  decaying  constituents.  Fluid  density  is  related  to  salinity,  in  practical  salinity  units 
(psu),  using  the  following  equation  of  state: 


p  =  p/+  %c 

where  pyis  the  reference  fluid  density  (that  is,  the  density  of  freshwater)  [ML-3],  and  3p /dC  is  the  slope  of  a  linear  relation  between 
fluid  density  and  salinity  [L0].  For  salinities  ranging  between  freshwater  and  typical  seawater,  dp/dC  has  an  approximate  value 
of  0.7.  The  effect  of  temperature  on  fluid  density  is  not  considered  here,  although  it  could  be  important  for  some  applications.  For 
the  Everglades  application,  seasonal  temperature  variations  can  be  substantial,  but  spatial  variations  are  assumed  to  have  a  negligi¬ 
ble  effect  on  flow.  Simultaneous  solutions  to  equations  1  through  5  result  in  spatial  distributions  for  h,  C,  p,  vx,  and 

2.2  Ground-Water  Flow  and  Solute  Transport 

Simulation  of  ground- water  flow  in  an  aquifer  with  spatially  varying  fluid  density  requires  solving  the  three-dimensional,  cou¬ 
pled  ground- water  flow  and  solute-transport  equations.  The  assumption  of  shallow  depths  (used  for  surface-water  flow)  does  not 
apply  to  ground  water,  and  a  full,  three-dimensional  solution  is  required  to  account  for  vertical  variations  in  aquifer  properties  and 
flow  patterns.  Guo  and  Langevin  (2002)  derive  a  variable-density  form  of  the  fully  saturated,  three-dimensional  ground- water  flow 
equation  in  terms  of  fy,  which  is  equivalent  freshwater  head  [L]  (Lusczynski,  1961): 


dp  dC 
dC  dt 


~Psg<lsg+Pflr+Pflet 


(6) 


where  K^xx,  K^yy  and  Kfzz  are  equivalent  freshwater  hydraulic  conductivities  [LT1]  in  the  x-9  y-,  and  z-directions,  psg  is  the  density 
of  the  leakage  fluid  calculated  by  substituting  Csg  into  equation  5  [ML-3],  Sf  is  the  specific  storage  in  terms  of  equivalent  freshwater 
head  [L_1],  and  0  is  porosity  [L0].  The  governing  equation  for  solute  transport  within  a  porous  medium  (Zheng  and  Wang,  1999) 
is  written  as: 
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where  qx,  q},  and  qz  are  the  specific  discharges  in  the  x-,  >>-,  and  "-directions  [LT1].  Equations  6  and  7  are  coupled  in  two  ways. 
First,  the  fluid  density  terms  in  equation  6  are  related  to  solute  concentrations  through  the  equation  of  state  (eq.  5).  Second,  the 
solute-transport  equation  (eq.  7)  contains  specific  discharge  terms  ( qx ,  qy  and  qz)  that  result  from  a  solution  to  the  ground- water 
flow  equation  (eq.  6). 
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2.3  Surface-Water/Ground-Water  Interactions 

A  variety  of  methods,  and  combinations  thereof,  were  evaluated  for  calculating  flow  interactions  between  surface  water  and 
ground  water,  including  a  simple  Darcy  equation,  modified  versions  of  the  Green- Ampt  infiltration  equation  (Green  and  Ampt, 
1911),  and  a  solution  to  the  Richard’s  equation  (Richards,  1931).  Field  observations  and  model  results  confirm  that  unsaturated 
zones  are  rarely  encountered  in  the  Everglades  coastal  wetlands,  but  if  encountered,  they  are  thin  and  of  short  duration.  Based  on 
these  observations,  a  simple  variable-density  form  of  Darcy’s  equation  (Juster,  1995;  Guo  and  Langevin,  2002)  written  in  terms  of 
equivalent  freshwater  head  was  programmed  to  calculate  vertical  leakage  between  the  wetlands  and  aquifer.  If  a  thin  unsaturated 
zone  develops  during  the  simulation,  leakage  rates  are  constrained  such  that  rates  do  not  increase  as  the  water  table  drops  farther 
below  land  surface  (described  in  the  next  section).  The  leakage  flux  is  applied  as  a  source/sink  term  in  the  continuity  equation  for 
surface-water  flow  (eq.  1)  and  as  a  boundary  flux  to  the  aquifer  surface  for  the  ground- water  system.  The  difference  in  treatment 
is  due  to  a  two-dimensional  surface-water  model  and  a  three-dimensional  ground- water  model.  Vertical  leakage  is  calculated  using 
the  following  variable-density  form  of  Darcy’s  law: 
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In  this  paper,  solute  mass  transfer  between  systems  is  assumed  to  occur  solely  through  advection.  Thus,  the  solute  mass  flux  is 
simply  the  product  of  leakage  and  the  upstream  solute  concentration  of  the  leakage  fluid. 


3.  Numerical  Implementation 

To  solve  the  coupled  surface-water  and  ground- water  equations  presented  in  the  previous  section,  the  finite-difference  pro¬ 
grams,  SWIFT2D  and  SEAWAT,  were  modified  to  run  timesteps  sequentially  under  the  control  of  a  master  program  called 
FTLOADDS  (Flow  and  Transport  in  a  Linked  Overland/Aquifer  Density  Dependent  System).  The  SWIFT2D  two-dimensional 
hydrodynamic  flow  and  solute-transport  code  was  originally  developed  for  bays  and  shallow  estuaries  (Leendertse  and  Gritton, 
1971;  Leendertse  1987).  The  code  has  been  applied  to  Jamaica  Bay,  New  York  (Leendertse,  1972),  to  Delta  Works,  The  Nether¬ 
lands  (Leendertse  and  others,  1981),  to  Tampa,  Florida  (Goodwin,  1987;  Goodwin  1991),  to  Pamlico  River  Estuary,  North  Carolina 
(Bales  and  Robbins,  1995),  to  Charlotte  Harbor,  Florida  (Goodwin,  1996)  and  to  the  Neuse  River  Estuary,  North  Carolina  (Robbins 
and  Bales,  1995).  The  SWIFT2D  program  was  later  modified  by  Swain  and  others  (2004)  to  represent  overland  flow  in  coastal 
wetlands  and  to  include  the  effects  of  spatially  distributed  rainfall  and  evapotranspiration.  SWIFT2D  uses  a  finite-difference 
approximation  to  solve  the  surface-water  equations  (eqs.  1-5).  SEAWAT,  a  combined  version  of  MODFLOW  (McDonald  and  Har- 
baugh,  1988)  and  MT3DMS  (Zheng  and  Wang,  1999),  was  designed  to  solve  the  three-dimensional  variable-density  ground- water 
flow  and  solute-transport  equations  (eqs.  5-7)  using  finite-difference  methods  (Guo  and  Bennett,  1998;  Guo  and  Langevin,  2002; 
Langevin  and  others,  2003).  Examples  of  SEAWAT  applications  include  simulation  of  submarine  ground-water  discharge  to  a 
marine  estuary  (Langevin,  2001,  2003)  and  intercode  comparisons  (Bakker,  2003;  Bakker  and  others,  2004). 

SWIFT2D  uses  an  alternating-direction  implicit  (ADI)  method  and  a  space-  and  time- staggered  grid  to  solve  the  governing 
equations,  such  that  each  surface-water  timestep  is  divided  into  two  half  timesteps — one  half  timestep  for  flow  and  transport  in  the 
v-direction  and  the  other  for  the  y-direction.  In  each  of  the  two  phases  of  the  ADI  method,  the  continuity  equation  and  one  of  the 
components  of  the  momentum  equations  are  solved  with  local  storage  (and  corresponding  transport  term  of  the  continuity  equa¬ 
tion),  local  acceleration,  pressure  gradient,  and  the  frictional  term  of  the  momentum  equation  treated  implicitly.  The  last  three  terms 
on  the  left-hand  side  of  equation  1  (the  source  and  sink  terms)  are  not  included  in  the  finite-difference  solution,  but  are  separately 
added  to,  or  subtracted  from,  the  cell  volume.  SEAWAT  uses  an  implicit  finite-difference  approximation  to  solve  the  ground- water 
flow  equation  (eq.  6),  and  contains  several  alternative  methods  for  solving  the  solute-transport  equation  (eq.  7),  including  implicit 
and  explicit  finite-difference  methods  with  various  weighting  options  and  the  method  of  characteristics. 

The  integrated  code  for  SWIFT2D  and  SEAWAT  requires  cells  that  coincide  and  are  identical  in  size.  The  integrated  code 
was  designed  such  that  the  domains  of  the  two  models  need  not  be  identical,  provided  that  leakage  fluxes  can  be  neglected  in  areas 
where  the  two  models  do  not  overlap.  Although  not  used  for  the  Everglades  application,  this  feature  may  prove  useful  for  certain 
applications  where  the  extension  of  the  model  domain  is  necessary  in  only  one  of  the  two  systems. 
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3.1  Coupling  Procedure 


Panday  and  Huyakorn  (2004)  discuss  several  options  for  coupling  surface  and  subsurface  models:  (1)  a  “fully  coupled”  or 
“fully  implicit”  approach,  (2)  a  sequentially  coupled  approach  in  which  the  interaction  flux  is  applied  as  a  boundary  condition  to 
each  model,  or  (3)  a  sequentially  coupled  approach  in  which  the  head  for  one  system  acts  as  a  general-head  boundary  for  the  other 
system.  Fairbanks  and  others  (2001)  demonstrate  that  the  fully  implicit  approach,  in  which  a  single  set  of  matrix  equations  is  for¬ 
mulated  for  both  systems,  is  the  most  robust  and  accurate.  Reported  applications  with  the  “fully  coupled”  approach  have  been  lim¬ 
ited  to  using  a  diffusion  analogy  or  kinematic  wave  approximation  for  the  overland  flow  system.  The  sequentially  coupled 
approaches  may  be  programmed  to  use  an  iterative  coupling  scheme,  in  which  solutions  are  repeated  for  the  same  timestep  until 
the  change  between  subsequent  interaction  fluxes  is  less  than  a  user- specified  value,  or  a  time-lagged  approach.  One  advantage  of 
the  sequentially  coupled  approach  (used  here)  is  that  many  sub-timesteps  can  be  used  for  the  surface  model  before  solving  for  a 
longer  subsurface  timestep  (Fairbanks  and  others,  2001).  This  advantage  is  particularly  useful  for  the  Everglades  application, 
where  surface-water  timesteps  are  constrained  far  more  severely  than  ground- water  timesteps  because  of  rapid  surface-water  wave 
propagation  speeds. 

A  sequentially  coupled  time-lagged  approach  was  implemented  to  couple  the  surface-water  and  ground- water  systems.  The 
approach  is  mass  conservative  in  that  the  exact  leakage  flux  imposed  on  the  surface-water  system  is  also  imposed  on  the  ground- 
water  system.  Figure  2  shows  a  schematic  of  the  coupling  approach  for  a  single  ground- water  stress  period  (m),  which  is  a  time 
period  when  hydraulic  stresses  are  assumed  constant.  First,  surface-water  flow  is  simulated  for  each  sub-timestep  (n)  by  applying 
a  leakage  rate  calculated  with  the  ground- water  head  from  the  end  of  the  previous  stress  period  and  the  surface-water  level  at  the 
current  sub-timestep.  To  ensure  conservation  of  fluid  mass,  individual  leakage  quantities  for  each  surface  sub-timestep  are  summed 

according  to  the  following  equation  to  calculate  a  time- weighted  average  leakage  rate,  q™g  ,  for  the  stress  period: 


nsub 

I K4. 


m,n 

sg 


Vsg  = 


nsub 

X  A  t 

n  -  1 


(9) 


where  m  is  stress  period  number,  nsub  is  the  number  of  sub-timesteps  in  the  stress  period,  and  At  is  the  sub-timestep  length.  This 
average  leakage  rate  is  then  applied  to  the  ground-water  model  as  a  specified-flux  boundary  that  remains  constant  for  the  stress 
period. 


3 2  Leakage  Calculation 

The  leakage  flux  is  calculated  in  one  of  three  ways,  depending  on  the  presence  of  standing  surface  water  and  the  vertical 
position  of  the  water  table.  For  sub-timesteps  with  a  dry  surface-water  cell,  the  leakage  flux  is  set  to  zero: 

For  dry  surface-water  conditions, 
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Simulate  surface-water  flow  and  transport  using 
leakage  rate  calculated  with  surface-water  head 
at  current  sub-timestep  and  ground-water  head 
from  end  of  previous  stress  period 


Calculate  average  leakage  rate  for  stress  period 
by  summing  total  leakage  volume  over  all  sub- 
timesteps  and  dividing  by  stress  period  length 


Simulate  ground-water  flow  and  transport  using 
average  leakage  rate  distributed  evenly  over 
stress  period 


Surface-water 

sub-timesteps 


Time-weighted  average 
^  sg  leakage  rate 


stress  period  m 


Figure  2.  Relation  between  surface-water  and  ground-water  timesteps. 
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For  cells  with  standing  surface  water,  the  leakage  formulation  is  based  on  a  variable-density  form  of  Darcy’s  law,  where  the  con¬ 
ductance  term  is  calculated  using  the  mean  hydraulic  conductivity  between  land  surface  and  the  vertical  center  of  the  top  aquifer 
layer  (fig.  3).  Thus,  the  leakage  flow  length  is  from  land  surface  to  the  center  of  the  top  aquifer  layer,  where  the  equivalent  fresh¬ 
water  head  is  calculated  by  SEAWAT.  The  formulation  allows  for  the  presence  of  a  thin,  hydraulically  resistive  layer  at  the  land 
surface,  which  for  the  Everglades  application,  corresponds  to  the  peat  and  marl  unit  overlying  the  Biscayne  aquifer.  Thus,  for  cells 
where  the  surface-water  level,  hfjQ,  is  above  land  surface,  and  the  water-table  elevation,  is  above  the  bottom  of  the  thin 
layer  (as  shown  in  figure  3),  the  leakage  flux  is  calculated  using  the  following  variable-density  form  of  Darcy’s  law: 
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where  ZqTL  bot  *s  the  elevation  of  the  thin  layer  bottom  [L],  KfQj/2  is  the  thickness- weighted  harmonic  mean  average  of  equiv¬ 
alent  freshwater  hydraulic  conductivity  between  land  surface  and  the  center  of  aquifer  layer  1  [LT1],  AZ- ^  j  is  the  layer  1  cell  thick¬ 
ness,  hf  j  f  o  is  the  equivalent  freshwater  head  [L]  of  the  surface  water  (layer  0)  evaluated  at  land  surface,  hf  ^  i  is  the  equivalent 
freshwater  head  [L]  at  the  vertical  center  of  layer  1,  Pijj/2  is  the  average  surface-  and  ground-water  fluid  density  [ML-3],  and 
is  the  center  elevation  [L]  for  layer  1 . 


Figure  3.  Diagram  of  a 
surface-water  cell  and 
an  uppermost  aquifer 
cell.  Piezometers  are 
used  to  demonstrate  the 
concept  of  equivalent 
freshwater  head  and  the 
reference  elevations  for 
calculating  equivalent 
freshwater  head. 
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EXPLANATION 

hu1  is  the  head  in  a  piezometer  open  at  the  center  of  aquifer  layer  1 

hflJ1  is  the  equivalent  freshwater  head  in  a  piezometer  open  at  the  center  of  aquifer  layer  1 

hIJ0  is  the  elevation  of  the  water  surface 

hflj0  is  the  equivalent  freshwater  head  in  a  piezometer  open  at  land  surface 
CiJ0  is  the  solute  concentration  of  the  surface  water 
pIJ0  is  the  density  of  the  surface  water 

C:j1  is  the  solute  concentration  of  the  ground  water  in  aquifer  layer  1 

pu1  is  the  density  of  the  ground  water  in  aquifer  layer  1 

ZIJ0  is  the  elevation  of  land  surface 

ZI]Tl_bot  is  the  elevation  at  the  base  of  the  thin  layer 

AZ,  ,,  is  the  thickness  of  the  thin  layer 

Z,Jf1  is  the  center  elevation  of  aquifer  layer  1 

AZ;  ,  is  the  thickness  of  aquifer  layer  1 
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Although  rare  for  the  Everglades  application,  a  thin  unsaturated  zone  can  develop  if  the  surface-water  layer  is  rapidly  flooded 
while  the  water  table  remains  below  the  bottom  of  the  thin  hydraulically  resistive  layer.  If  these  conditions  occur  in  the  model, 
equation  10  is  no  longer  valid,  and  it  is  assumed  that  the  entire  head  loss  between  land  surface  and  the  center  of  the  cell  occurs 
across  the  thin  layer.  With  this  assumption  and  the  assumption  that  the  pressure  at  the  bottom  of  the  thin  layer  is  atmospheric,  the 
following  equation  is  used  to  approximate  the  flux  through  the  thin  unsaturated  zone: 

For  hi,j,0  >  Zi,j,0  and  hi,j,l  <  Zi,j,TL_BOT’ 
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where  K^i  j  TL  is  the  vertical  equivalent  freshwater  hydraulic  conductivity  of  the  thin  layer  [LT1],  A Zt  j  TL  is  the  thickness  of  the  thin 
layer  [L],  and  p  ^  is  the  density  of  the  surface  water  [ML-3].  This  approximation  for  the  flux  through  the  unsaturated  zone  is  based 
on  the  approach  used  by  MODFLOW  (McDonald  and  Harbaugh,  1988)  and  SEAWAT  (Guo  and  Langevin,  2002)  for  the  River 
package.  A  more  sophisticated  approach,  such  as  one  based  on  a  modified  Green-  Ampt  formulation  may  be  required  if  future  appli¬ 
cations  result  in  ponded  surface  water  overlying  a  relatively  thick  unsaturated  zone. 


3.3  Surface/Subsurface  Solute  Exchange 

A  mass-conservative  approach  was  designed  to  allow  advective  transport  between  the  wetland  and  underlying  aquifer.  Calcu¬ 
lation  of  the  advective  mass  flux  is  straightforward  for  most  stress  periods  in  which  the  advective  flux  is  either  up  or  down  for  the 
entire  period.  If  leakage  is  downward,  the  advective  mass  flux  is  the  product  of  the  leakage  rate  and  the  surface-water  solute  con¬ 
centration.  Likewise,  for  upward  leakage,  the  advective  mass  flux  is  the  product  of  the  leakage  rate  and  the  ground- water  solute 
concentration. 

Concentration  changes  that  result  from  advective  transport  between  the  wetland  and  aquifer  are  calculated  for  each  sub- 
timestep  in  SWIFT2D.  To  account  for  advective  leakage  transport  in  the  ground-water  model,  the  total  solute  mass  added  to,  and 
subtracted  from,  the  surface-water  cell  is  summed  in  SWIFT2D  for  each  stress  period.  This  total  mass  transferred  is  then  divided 
by  the  total  leakage  volume  for  that  stress  period  to  calculate  an  effective  leakage  concentration.  Thus,  the  volumetric  leakage  rate 
applied  to  SEAWAT  has  an  associated  effective  concentration  that  results  in  the  conservation  of  mass  between  the  two  systems.  If 
multiple  reversals  of  leakage  direction  occur  during  a  single  stress  period,  the  effective  concentration  can  be  very  small  or  even 
negative.  If  a  highly  saline  ground- water  system  is  overlain  by  a  fresh  surface-water  system  and  leakage  reverses  direction  multiple 
times,  then  the  net  salt  transfer  will  be  upward  even  if  the  net  leakage  is  downward.  In  this  case,  the  effective  concentration  will 
be  negative  indicating  that  the  net  salt  flux  is  in  the  opposite  direction  of  the  net  volume  flux. 


3.4  Spatially  Distributed  Rainfall  and  Evapotranspiration 

The  original  SWIFT2D  program  was  modified  to  include  spatially  distributed  rainfall  and  evapotranspiration  (Swain  and  oth¬ 
ers,  2004).  For  conditions  with  standing  surface  water,  rainfall  is  applied  to  the  surface-water  layer  with  a  solute  concentration  of 
zero.  If  a  surface-water  cell  is  dry,  the  rainfall  volume  is  applied  directly  and  instantaneously  to  the  water  table  in  layer  1  of  the 
ground-water  model.  The  same  approach  is  used  to  determine  where  the  evapotranspiration  flux  is  applied.  In  the  current  Ever¬ 
glades  application,  however,  the  evapotranspiration  flux  is  calculated  by  the  model  during  the  simulation  using  a  modified  Priestly- 
Taylor  approximation  to  the  physics-based  Penman-Monteith  model  as  described  by  German  (2000a;  2000b)  and  Swain  and  others 
(2004).  Evapotranspiration  rates  are  calculated  as  a  function  of  solar  radiation  and  water  depth.  Two  sets  of  coefficients  were  esti¬ 
mated  through  linear  regression  and  used  in  the  Priestly-Taylor  approximation.  One  set  was  calculated  for  conditions  with  standing 
water;  the  other  was  calculated  for  dry  surface-water  conditions  when  the  water  table  was  below  land  surface.  For  conditions  with 
standing  water,  the  entire  evapotranspiration  flux  is  withdrawn  from  surface  water,  rather  than  withdrawing  evaporation  from  sur¬ 
face  water  and  transpiration  from  ground  water. 

3.5  Rewetting 

In  the  present  integrated  code,  surface-water  and  ground- water  cells  are  allowed  to  dry  and  rewet.  The  original  versions  of 
SWIFT2D  and  SEAWAT  both  have  options  for  cells  to  rewet  from  the  four  surrounding  nodes,  and  in  the  case  of  SEAWAT,  from 
an  underlying  node.  Rewetting  in  SEAWAT  is  a  simple  extension  of  the  procedure  implemented  in  MODFLOW  (McDonald  and 
others,  1992).  In  the  integrated  program,  a  modification  was  made  to  SWIFT2D  to  allow  surface-water  cells  to  rewet  from  the 
underlying  ground-water  cell  if  the  water  table  rises  above  land  surface.  If  a  surface-water  cell  is  dry,  then  SWIFT2D  compares 
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surrounding  stages  and  height  of  the  water  table  with  land-surface  elevation.  This  comparison  is  performed  for  a  sub-timestep  inter¬ 
val  provided  by  the  user.  If  a  surrounding  stage  or  water-table  elevation  is  above  land  surface,  the  cell  is  reactivated  and  included 
in  the  computational  domain  for  the  subsequent  sub-timestep. 


4.  Application  of  Integrated  Model  to  the  Southern  Florida  Everglades 

The  integrated  model  was  applied  to  the  southern  Everglades  and  northeastern  Florida  Bay  (fig.  4)  to  evaluate  the  dominant 
hydrologic  processes,  including  surface-water  and  ground- water  interactions,  and  to  synthesize  a  wide  range  of  hydrologic  data 
collected  for  the  area.  The  specific  objective  of  the  model  application  was  to  develop  a  numerical  tool  that  could  be  used  to  quantify 
freshwater  discharges  to  northeastern  Florida  Bay,  predict  temporal  and  spatial  variations  in  coastal  salinity  patterns,  and  represent 
wetland  hydroperiods.  Presently,  the  numerical  tool  is  being  used  to  evaluate  the  effects  of  the  Comprehensive  Everglades  Resto¬ 
ration  Plan  (CERP)  on  future  hydrologic  conditions  (heads,  flows,  and  salinities)  in  the  coastal  wetlands  and  adjacent  Florida  Bay 
estuary.  The  principal  purpose  of  CERP  is  to  restore  the  southern  Florida  ecosystem,  which  includes  the  Everglades  (http:// 
w  w  w .  e  vergladesplan .  org) . 
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Figure  4.  Study  area  including  location  of  Taylor  Slough,  Florida  Bay,  monitoring  stations,  and  model  boundary. 
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The  U.S.  Geological  Survey  began  modeling  Taylor  Slough  hydrology  in  1995  with  the  development  of  a  hydrodynamic  sur¬ 
face-water  flow  and  solute-transport  model,  as  described  by  Swain  and  others  (2004).  Information  relevant  to  the  application  of 
the  present  integrated  code  is  briefly  summarized  here.  The  surface-water  model  by  Swain  and  others  (2004)  simulated  the  1-year 
period  from  August  1996  to  July  1997.  As  treatment  of  surface-water  and  ground-water  exchanges  was  not  included  in  the  early 
model,  the  lack  of  a  ground- water  component  was  a  limitation,  considering  the  documented  importance  of  surface-water  and 
ground-water  exchanges  in  southern  Florida  (for  example,  Merritt,  1996a;  Swain  and  others,  1996).  Subsequently,  the 
FTLOADDS  program  was  designed  and  a  SEAWAT  ground-water  model  was  developed  for  the  study  area.  The  integrated  sim¬ 
ulation,  described  herein,  was  extended  to  7  years  to  represent  the  period  from  January  1996  to  December  2002.  There  are  a  total 
of  2,542  one-day  ground- water  stress  periods,  with  each  day  divided  into  7.5-minute  surface-water  sub-timesteps.  The  model  also 
includes  a  15-day  “warm-up  period”  in  which  only  the  surface-water  system  is  simulated.  The  initial  condition  for  the  simulation 
is  a  fresh,  flat  pool  at  an  elevation  of  0.78  m  above  the  North  American  Vertical  Datum  of  1988  (NAVD  88).  Leakage  and  solute 
mass  transfer  to  or  from  the  surface-water  system  is  included  during  the  warm-up  period  using  initial  ground-water  heads  and  salin¬ 
ities  that  remain  constant  during  the  15 -day  period.  Initial  ground- water  heads  and  salinities  were  set  using  the  results  from  a  pre¬ 
liminary  7-year  simulation  that  had  been  rerun  seven  times  until  the  model  reached  dynamic  equilibrium. 


4.1  Description  of  Study  Area 

The  focus  of  the  model  application  is  a  900-km2  area  of  southeastern  Everglades  National  Park  that  includes  northeastern  Flor- 
ida  Bay  (fig.  4).  The  study  area  encompasses  Taylor  Slough,  which  is  the  smaller  of  two  main  Everglades  sloughs  in  southern  Flor¬ 
ida.  Rainfall  is  a  dominant  source  of  freshwater  for  the  study  area,  which  receives  an  average  of  about  140  cm/yr.  Inflows  also 
occur  by  means  of  a  water-management  system  that  controls  water  levels  in  southern  Florida  to  prevent  flooding.  A  pump  structure 
(S332;  not  shown)  is  used  to  transfer  water  from  the  Levee  3 1 W  (L-3 1 W)  canal  into  the  wetlands  just  north  of  Taylor  Slough  Bridge 
(fig.  4).  A  portion  of  this  water  then  flows  beneath  Taylor  Slough  Bridge  into  the  main  part  of  Taylor  Slough,  which  extends  south 
to  Little  Madeira  Bay  and  Florida  Bay  and  east  toward  Joe  Bay.  The  L-31W  canal,  which  extends  southward  from  the  northern 
model  boundary,  is  another  contributor  of  freshwater  to  the  coastal  wetlands  and  usually  flows  only  during  flood  control  operations. 
Inflows  are  also  common  from  the  C-l  1 1  canal  between  control  structures  S18C  and  SI 97.  A  continuous  levee  once  existed  on  the 
south  side  of  C-l  1 1  between  S18C  and  SI 97.  In  the  1970’ s,  notches  were  cut  through  the  levee  to  allow  surface  water  from  C-l  1 1 
to  flow  into  the  coastal  wetlands.  In  1998,  much  of  the  remaining  levee  was  degraded  to  increase  inflows  into  the  wetland  area 
south  of  S18C.  On  the  west  side  of  the  study  area,  culverts  beneath  the  Main  Park  Road  allow  surface  water  to  exchange  with  the 
wetlands  to  the  west,  but  flow  measurements  indicate  the  exchanges  are  minimal  (Tillis,  2001;  Stewart  and  others,  2002) 

A  principal  hydrologic  feature  in  the  study  area  is  the  Buttonwood  Embankment — a  nearly  continuous  ridge  along  the  Florida 
Bay  coastline.  This  ridge  is  observed  to  be  about  0.3  m  higher  than  the  surrounding  marsh,  and  was  formed  either  by  the  buildup 
of  organic  detritus  from  the  stands  of  mangrove  forest  or  by  sediment  deposition  from  Florida  Bay  during  periodic  hurricanes  and 
tropical  storms  (Holmes  and  others,  2000).  The  ridge  itself  forms  a  partial  low-crowned  barrier  and  obstructs  direct  overland  flow 
from  the  coastal  wetlands  into  northeastern  Florida  Bay  at  most  times.  Hydraulic  connection  between  the  coastal  wetlands  and 
northeastern  Florida  Bay  occurs  through  coastal  creeks  that  have  incised  the  Buttonwood  Embankment.  Overtopping  is  infrequent 
and  is  typically  caused  by  northward  moving  storms  or  hurricanes  that  force  brackish  Florida  Bay  water  over  the  embankment  and 
into  the  coastal  wetlands  (Hittle,  2000). 

In  1995,  the  U.S.  Geological  Survey  began  collecting  continuous  (15-minute  interval)  discharge,  stage,  and  salinity  data  at  the 
mouths  of  coastal  creeks  (Hittle  and  others,  2001).  These  field  data  quantify  exchange  rates  between  the  coastal  wetlands  and  north¬ 
eastern  Florida  Bay  and  provide  insight  about  the  hydrologic  processes  driving  the  exchanges.  The  field  data  reveal  the  following 
general  patterns  and  trends  (Hittle,  2000):  (1)  the  average  annual  net  discharge  to  Florida  Bay  through  the  five  measured  coastal 
creeks  (McCormick,  Mud,  Trout,  Taylor  River,  and  West  Highway)  is  3.2  x  108  m3/yr,  with  about  80  percent  of  the  discharge 
occurring  during  the  summer  months  from  May  to  October;  (2)  most  of  the  discharge  (60  percent)  to  northeastern  Florida  Bay 
occurs  through  Trout  Creek;  (3)  wind  variations  seem  to  cause  flow  reversals,  forcing  brackish  Florida  Bay  waters  up  into  the  bays 
and  coastal  wetlands;  and  (4)  the  interconnected  mud  banks  of  Florida  Bay  and  the  Florida  Keys  substantially  dampen  the  tidal 
signature  in  the  northeastern  part  of  the  bay. 

Field  studies  of  surface-water  and  ground- water  interactions  within  the  Taylor  Slough  area  are  reported  by  Harvey  and  others 
(2000a,  2000b)  and  Price  (2001).  Using  a  chloride  dilution  method,  Harvey  and  others  (2000a)  indicated  that:  (1)  upward  ground- 
water  flow  in  November  1997  may  have  been  as  large  as  3  cm/d  in  the  area  near  NP67  and  TSH  (fig.  4),  (2)  the  high  water  levels 
on  the  northwestern  side  of  Old  Ingraham  Highway  are  probably  driving  ground- water  flow  into  the  western  part  of  Taylor  Slough, 
and  (3)  upward  ground- water  flow  within  the  slough  itself  is  about  0.06  cm/d  for  the  period  between  September  1997  and  Septem¬ 
ber  1999.  Using  a  variety  of  geochemical  tracers,  Price  (2001)  estimated  leakage  rates  that  ranged  over  4  orders  of  magnitude. 

Salinity  patterns  in  northeastern  Florida  Bay  are  highly  influenced  by  freshwater  discharges  from  the  coastal  wetlands  (Hittle, 
2000;  Nuttle  and  others,  2000).  Flows  from  coastal  creeks  appear  to  be  the  dominant  source  of  freshwater  for  northeastern  Florida 
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Bay,  although  some  scientists  have  suggested  that  overtopping  of  the  Buttonwood  Embankment  and  ground  water  are  important 
contributors.  Substantial  freshwater  inputs  into  Florida  Bay  from  shallow  ground  water  seem  unlikely  because  recent  helicopter 
electromagnetic  surveys  revealed  that  shallow  saline  ground  water  extends  at  least  7  km  inland  (fig.  4;  Fitterman  and  Deszcz-Pan, 
1998).  Corbett  and  others  (1999)  used  geochemical  tracers  to  identify  submarine  ground- water  discharge  patterns  in  Florida  Bay 
and  discovered  what  appeared  to  be  higher  discharge  values  near  the  Florida  Keys  than  just  south  of  the  Buttonwood  Embankment. 
They  also  reported  that  a  seepage  meter  placed  in  northern  Florida  Bay  recorded  leakage  rates  in  excess  of  1  cm/d. 


4.2  Model  Design 

Application  of  the  integrated  model  to  the  southern  Everglades  and  northeastern  Florida  Bay  required  a  wide  variety  of  input 
data  for  both  the  surface-water  and  ground-water  systems.  Fortunately,  the  study  area  has  been  the  focus  of  concentrated  hydro- 
logic  investigations.  Many  of  these  investigations  were  initiated  primarily  to  collect  data  for  the  integrated  model.  Brief  descrip¬ 
tions  are  provided  here  for  many  of  the  important  hydrologic  input  parameters,  such  as  land-surface  elevation,  evapotranspiration, 
and  aquifer  hydraulic  conductivity. 

The  finite-difference  model  grid  used  for  the  Everglades  application  consists  of  98  east- west  rows  and  148  north- south  col¬ 
umns.  Model  cells  are  square  and  measure  305  m  per  side.  The  area  of  the  model  that  corresponds  to  the  coastal  wetlands  (the 
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active  model  area  north  of  the  Buttonwood  Embankment)  is  about  6.2  x  10  m  .  Florida  Bay  comprises  about  2.8  x  10  m  of 
the  model  area.  Fand- surface  elevations  were  calculated  for  each  cell  (fig.  5)  by  interpolating  values  from  a  helicopter  global 
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Figure  5.  Grid  of  land-surface  elevation,  relative  to  the  North  American  Vertical  Datum  of  1988  (NAVD  88),  as 
used  in  the  numerical  model. 
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positioning  system  (GPS)  survey  with  400-m  spacing  between  measurements  (Desmond,  2003).  Data  from  National  Oceanic  and 
Atmospheric  Administration  (NO  A  A)  nautical  charts  and  from  Hansen  and  Dewitt  (2000)  were  used  to  assign  bathymetry  for  Flor¬ 
ida  Bay  and  its  subembayments.  The  Buttonwood  Embankment  was  included  in  the  model  using  the  barrier  feature  in  SWIFT2D 
with  a  specified  sill  height  of  0.3  m  above  land  surface.  The  three-dimensional  ground- water  model  has  10  layers.  The  top  layer 
extends  from  land  surface  to  an  elevation  of  3.2  m  below  NAVD  88.  The  remaining  layers  are  uniform  in  volume  and  have  a  con¬ 
stant  thickness  of  3.2  m.  The  ground- water  model  extends  from  land  surface  to  the  base  of  the  Biscayne  aquifer,  as  defined  by  Fish 
and  Stewart  (1991)  and  revised  by  Fitterman  and  others  (1999).  The  lower  part  of  the  surficial  aquifer  system,  as  described  by 
Jarosewich  and  Wagner  (1985)  and  Fish  and  Stewart  (1991),  is  not  included  in  the  model.  Within  the  study  area,  the  base  of  the 
Biscayne  aquifer  thickens  to  the  east.  Ground- water  model  cells  are  inactive  if  the  vertical  center  of  the  model  cell  is  located 
beneath  the  base  of  the  Biscayne  aquifer.  This  approach  treats  the  base  of  the  Biscayne  aquifer  as  a  no-flow  boundary — an  approach 
commonly  used  in  southern  Florida  (Merritt,  1996a,  1996b;  Swain  and  others,  1996;  Langevin,  2001)  and  justified  by  the  sharp 
contrast  in  permeabilities. 

Boundary  conditions  for  the  surface-water  model  were  specified  for  the  model  perimeter  based  on  the  presence  of  roads, 
canals,  culverts,  islands,  and  an  estimated  sufficient  offshore  distance  from  the  southern  Florida  coastline  (fig.  6).  The  type  of 
boundary  used  for  each  segment  in  figure  6  is  listed  in  table  1 .  Boundary  levels  and  fluxes  were  assigned  based  on  field  data  using 
the  highest  temporal  frequency  available.  Two  different  types  of  boundaries  were  used  for  the  ground- water  model.  North  of  the 
Florida  Bay  coastline,  general-head  boundaries  were  applied  to  each  layer  of  the  ground- water  model  based  on  interpolated  head 
values  from  nearby  surface-water  and  ground- water  monitoring  sites  for  each  day  of  the  simulation  period.  Salinity  values  assigned 
to  the  general-head  boundaries  in  all  layers  were  estimated  from  the  airborne  geophysical  survey  (Fitterman  and  Deszcz-Pan,  1998). 
South  of  the  Florida  Bay  coastline  (corresponding  to  SW4,  SW5,  and  SW6  in  fig.  6),  a  no-flow  boundary  was  imposed  on  the 
ground- water  model.  A  general-head  boundary  was  originally  used  for  the  southern  ground- water  boundary,  but  field  data  were 
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Figure  6.  Model  grid  and  index  numbers  for  model  boundaries. 
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not  available  for  assignment  of  boundary  heads  and  concentrations.  Thus  a  no-flow  condition  was  assigned  based  on  the  assump¬ 
tion  that  the  boundary  is  far  south  (about  1  to  10  km)  of  the  area  where  ground- water  discharge  might  occur.  Simulated  leakage 
maps  confirmed  this  assumption. 

Table  1.  Description  of  boundary  conditions  for  the  surface-water  Rainfall  and  evapotranspiration  are  the  primary  sources  and 

model.  sinks  within  the  model  domain.  Rainfall  data  with  recorded 

[Locations  of  model  boundaries  are  in  figure  6.  D,  specified  discharge  bound-  intervals  as  short  as  15  minutes  were  spatially  interpolated 

ary;  NF,  no-flow  boundary;  WL,  specified  water-level  boundary]  (kriged)  for  each  model  cell  and  time  interpolated  for  each  sur- 

face- water  sub-timestep.  Evapotranspiration  was  included  in  the 
model  using  a  modified  form  of  the  Priestly  Taylor  approxima¬ 
tion  to  the  physics-based  Penman-Monteith  equation.  Site-spe¬ 
cific  evapotranspiration  data  from  a  study  by  German  (1999) 
were  used  to  calculate  coefficients  within  the  Priestly-Taylor 
approximation  such  that  evapotranspiration  rates  are  a  function 
of  water  depth  and  solar  radiation.  The  dependency  on  water 
depth  is  unusual,  but  appears  to  be  related  to  sheltering  by  vege¬ 
tation  and  submergence  of  vegetation  (German,  1999).  Separate 
coefficients  were  determined  for  periods  when  water  levels  were 
below  land  surface.  Swain  and  others  (2004)  provide  a  detailed 
description  for  the  methods  used  to  assign  rainfall  and  evapo¬ 
transpiration  to  the  model. 

The  surface-water  simulation  is  controlled  by  spatially 
varying  parameters  that  represent  relevant  processes.  This 
includes  defining  the  frictional  resistance  to  flow,  wind  friction 
factor  and  sheltering  coefficient,  and  dispersion  coefficient.  The 
frictional  resistance  to  flow  is  expressed  with  Manning’s  coeffi¬ 
cients.  Because  of  the  importance  of  this  term,  field  and  laboratory  research  was  performed  to  determine  the  effective  frictional 
resistance  to  water  flow  through  differing  Everglades  vegetation  types  (Lee  and  others,  1999).  Extensive  hydraulic  measurements 
(velocity,  depth,  hydraulic  gradient,  and  so  forth)  were  made  in  a  laboratory  flume  containing  transplanted  marsh  vegetation.  Field 
measurements  of  velocity,  depth,  and  vegetation  type  and  density  were  also  made  in  conjunction  with  point  measurements  of  the 
hydraulic  gradient  using  a  portable  pipe  manometer  at  many  locations  in  the  study  area  (Lee  and  others,  2000).  Results  of  these 
studies  indicate  high  Manning’s  n  values  and  relatively  small  variations  between  vegetation  types.  Values  of  Manning’s  n  vary  spa- 
daily  and  range  from  0.38  to  0.46  s/m  (Swain  and  others,  2004).  Open-water  areas  are  assigned  a  nominal  value  of  0.02  s/m  . 
Frictional  resistance  values  for  the  coastal  creeks  were  determined  by  calibration  and  from  field  measurements  at  two  monitoring 
stations  on  Taylor  River  (fig.  4).  The  measured  discharge  and  water-level  differences  between  the  original  station  and  the  upstream 
station  were  used  to  determine  n  from  Manning’s  equation.  The  computed  Manning’s  n  for  each  day  was  averaged  for  the  period 
and  yielded  a  mean  Manning’s  n  value  of  0.121  s/m1/3  (Swain  and  others,  2004).  The  calibrated  values  of  Manning’s  n  for  the 
coastal  creeks  ranged  from  0.058  to  0.152  s/m1/3.  The  coefficient  for  the  wind-friction  term  that  related  the  wind  velocity  squared 
to  the  rate  of  momentum  change  in  the  water  flow  has  a  specified  value  of  about  1.2  x  10"3  for  winds  less  than  36  m/s  (Large  and 
Pond,  1981).  This  coefficient  is  uniform  for  the  entire  study  area.  A  wind- sheltering  term  also  is  applied  to  account  for  the  effects 
of  emergent  vegetation.  Estimated  values  for  this  wind- sheltering  term  range  from  0.1  to  0.5  (Reid  and  Whitaker,  1976);  a  value 
of  0.33  is  used  in  the  model.  Jenter  and  Duff  (1999)  suggest  that  the  values  estimated  by  Reid  and  Whitaker  (1976)  are  reasonable 
for  the  Everglades  coastal  wetlands.  The  magnitude  of  the  dispersion  coefficient  for  surface  water  is  scale  dependent,  increasing 
with  the  size  of  the  water  body.  The  effective  dispersion  coefficient  is  on  the  order  of  1-10  m2/s  in  open  channels,  and  2  orders  of 
magnitude  greater  in  estuaries  (Fischer  and  others,  1979).  In  the  application  of  the  dispersion  coefficient  in  a  numerical  model,  the 
length  scale  of  importance  is  the  cell  size.  Therefore,  the  dispersion  coefficient  was  calibrated  by  matching  salinity  values  at  the 
coastal  creek  measurement  stations. 

Application  of  Darcy’s  law  to  calculate  leakage  rates  (eq.  9)  required  an  assignment  of  aquifer  properties  to  the  upper  half  of 
layer  1  of  the  ground- water  model.  The  current  program  reads  the  thickness  and  vertical  hydraulic  conductivity  of  the  thin  layer 
and  the  vertical  hydraulic  conductivity  of  the  part  of  the  aquifer  between  the  base  of  the  thin  layer  and  the  vertical  center  of  model 
layer  1 .  Peat  thickness  was  measured  at  74  locations  within  Taylor  Slough  by  pushing  a  steel  rod  through  the  peat  to  the  limestone 
surface  (Harvey  and  others,  2000b).  Reported  thickness  values  range  between  0.3  and  2.5  m  with  an  average  value  of  1.1  m.  Har¬ 
vey  and  others  (2000b)  also  measured  peat  hydraulic  conductivity  by  performing  bail  tests  at  seven  shallow  piezometer  sites. 
Reported  hydraulic  conductivity  values  from  that  study  are  0.09,  0.2,  0.2,  0.3,  0.4,  0.5,  and  1.3  m/d.  For  Florida  Bay,  leakage  coef¬ 
ficients  were  assigned  based  on  mapped  bottom  types  (Prager  and  Halley,  1997).  For  hard-bottom  areas  such  as  Joe  Bay,  a  vertical 
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aquifer  hydraulic  conductivity  value  of  0.75  m/d  was  used.  All  other  bottom  types  in  Florida  Bay  were  assumed  to  have  a  1-m  thick 
sediment  layer  with  a  vertical  hydraulic  conductivity  value  of  0. 1  m/d.  The  remaining  part  of  the  Biscayne  aquifer  was  assumed 
to  be  isotropic  and  homogeneous  with  vertical  and  horizontal  hydraulic  conductivity  values  of  0.75  and  5,000  m/d,  respectively. 
These  values,  which  were  determined  through  calibration,  compare  closely  with  values  used  in  other  numerical  models  of  the  area 
(Merritt,  1996a,  1996b;  Swain  and  others,  1996;  Langevin,  2001). 

Limited  data  exist  on  the  dispersive  properties  of  the  Biscayne  aquifer.  Merritt  (1996b)  developed  a  numerical  model  to  sim¬ 
ulate  a  chloride  plume  that  resulted  from  a  flowing  artesian  well  open  to  the  brackish  Floridan  aquifer.  The  simulated  plume  cov¬ 
ered  a  3  x  108  m2  area,  which  is  similar  in  scale  to  the  present  study  (9  x  108  m2).  Merritt  (1996b)  assigned  a  porosity  value  of  0.2 
based  on  specific  yield  estimates  for  the  porous  limestone  aquifer,  and  through  calibration,  estimated  values  of  76  and  0.03  m  for 
longitudinal  and  transverse  dispersivity,  respectively.  Based  on  the  relatively  low  value  for  transverse  dispersivity,  Merritt  (1996b) 
concluded  that  the  transverse  dispersion  simulated  by  the  model  was  caused  by  seasonal  variations  in  flow  velocity  and  direction, 
rather  than  by  mechanical  dispersion.  For  the  present  simulation,  hydrodynamic  dispersion  within  the  ground- water  model  was  not 
active  during  the  integrated  simulations.  Instead,  it  was  assumed  that  numerical  dispersion  resulting  from  the  solution  to  the  trans¬ 
port  equation  was  similar  in  magnitude  to  actual  dispersion.  This  assumption  was  supported  by  simulated  positions  of  the  fresh¬ 
water/saltwater  interface  in  the  Biscayne  aquifer  that  matched  observed  positions.  As  additional  data  on  the  interface  position, 
movement,  and  width  become  available,  the  model  may  be  refined  to  include  mechanical  dispersion. 


Table  2.  Calibration  statistics  for  daily  average  coastal  creek  dis¬ 
charges,  surface-water  stage  or  ground-water  head,  and  coastal 
creek  salinities. 

[Mean  errors  were  calculated  by  subtracting  measured  values  from  simulated 
values.  Station  locations  are  in  figure  4.  ME,  mean  error;  MAE,  mean  abso¬ 
lute  error;  RMSE,  root  mean  square  error;  Count,  number  of  data  points  used 
to  calculate  statistics;  m3/s,  cubic  meters  per  second;  psu,  practical  salinity 
units]. 


Station 

ME 

MAE 

RMSE 

Count 

Discharge  (m3/s) 

McCormick 

0.20 

1.56 

2.04 

2510 

Mud 

0.39 

1.86 

3.86 

2530 

Trout 

-1.78 

5.01 

7.06 

2526 

Taylor  River 

-0.23 

1.21 

2.90 

2554 

West  Highway 

-0.30 

1.10 

1.65 

2479 

Stage/Head  (m) 

NMP 

0.02 

0.02 

0.02 

2290 

CY3 

-0.07 

0.07 

0.07 

2275 

NP46 

0.00 

0.06 

0.09 

2475 

NP67 

0.01 

0.06 

0.08 

2493 

CY2 

-0.03 

0.04 

0.05 

2222 

TSH 

0.00 

0.06 

0.08 

2527 

E146 

0.05 

0.06 

0.09 

2477 

CHP 

-0.02 

0.06 

0.08 

2473 

EPSW 

0.08 

0.09 

0.10 

2461 

EVER  6 

-0.04 

0.06 

0.07 

2394 

EVER  7 

-0.04 

0.05 

0.07 

2444 

R127 

0.02 

0.07 

0.10 

2494 

P37 

0.00 

0.05 

0.07 

2465 

G-3619 

-0.04 

0.09 

0.12 

2438 

G-3353 

0.14 

0.15 

0.17 

2451 

G-1251 

0.04 

0.08 

0.10 

2034 

Salinity  (psu) 

McCormick 

2.76 

7.14 

9.43 

2508 

Mud 

2.10 

3.95 

5.08 

2421 

Trout 

2.33 

4.86 

6.45 

2529 

Taylor  River 

4.95 

6.35 

7.70 

2515 

West  Highway 

-1.43 

4.60 

5.57 

2512 

4.3  Model  Calibration  and  Analysis  of 
Simulation  Results 

Computer  runtimes  in  excess  of  30  hours  (on  a 
Pentium  IV  processor  running  at  1 .7  gigahertz)  for 
the  7-year  simulation  period  prohibited  use  of 
formalized  parameter  estimation  techniques. 
Instead,  calibration  was  achieved  by  judicious 
adjustment  of  the  input  parameters  that  seemed  to 
have  the  largest  uncertainty  range  and  the  largest 
effect  on  simulation  results.  Calibration  statistics 
for  coastal  creek  discharges,  wetland  stages,  and 
creek  salinities  are  reported  in  table  2.  In  addition  to 
the  model  results  given  in  table  2,  other  model 
results  were  also  compared  with  field  data.  For 
example,  the  simulated  aquifer  salinities  were 
carefully  evaluated  to  ensure  the  model  adequately 
matched  the  results  from  an  airborne 
electromagnetic  survey  (Fitterman  and  Deszcz-Pan, 
1998). 

A  water  budget  was  prepared  from  model 
results  for  the  coastal  wetland  part  of  the  model 
domain,  north  of  the  Florida  Bay  coastline  (fig.  7 
and  table  3).  The  water  budget  is  for  surface  water 
and  does  not  include  lateral  ground- water  flows  or 
evapotranspiration  directly  from  the  water  table, 
which  is  about  45  cm.  Water  budget  components  are 
given  as  annual  average  values  for  individual  years 
(table  3)  and  for  the  7-year  simulation  period  (fig.  7 
and  table  3).  Coastal  creek  discharges,  and  other 
discharge  values,  were  divided  by  the  wetland  area 
(6.2  x  108  m2)  to  give  length  units  that  can  be 
compared  directly  with  rainfall  and 
evapotranspiration. 
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JL  17.1  ANNUAL  AVERAGE  NET  INFLOW  TO  THE  COASTAL  WETLAND 
FOR  THE  7-YEAR  SIMULATION  PERIOD  (1996-2002)- 
Values  are  in  centimeters  per  year 


f  1 7.4  ANNUAL  AVERAGE  NET  OUTFLOW  FROM  THE  COASTAL 

WETLAND  FOR  THE  7-YEAR  SIMULATION  PERIOD  (1996-2002)- 
Values  are  in  centimeters  per  year 


Figure  7.  Annual  average  surface-water  budget  for  the  coastal  wetland  for  the  7-year  simulation  period  (1996-2002). 
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Table  3.  Simulated  surface-water  budget  for  the  coastal  wetland. 

[Values  are  expressed  in  centimeters  as  inflows  or  outflows  to  the  wetland.  Discrepancies  between  total  inflows  and  outflows  are  due  to 
changes  in  storage  and  model  error.] 


Inflows 

Name 

Year 

1996  1997  1998  1999  2000  2001  2002 

Average 

Annual 

Rainfall 

128.32 

159.21 

127.93 

167.63 

102.61 

143.25 

134.05 

137.57 

Leakage 

15.71 

16.18 

15.82 

17.48 

19.96 

14.05 

20.76 

17.14 

C-lll 

21.44 

26.93 

27.34 

24.47 

30.88 

26.21 

22.84 

25.73 

TSB 

9.74 

15.31 

11.74 

15.84 

15.87 

15.57 

11.56 

13.66 

L-31W 

3.06 

4.00 

3.26 

15.91 

3.72 

0.00 

0.28 

4.32 

Trout 

8.75 

9.28 

8.04 

10.27 

9.54 

9.18 

10.31 

9.34 

Mud 

1.24 

1.34 

1.69 

4.04 

2.99 

2.30 

2.62 

2.32 

West  Highway 

0.26 

0.32 

0.11 

0.57 

0.09 

0.41 

0.38 

0.30 

East  Creek 

0.97 

1.14 

1.10 

1.61 

1.34 

1.55 

1.21 

1.27 

Taylor  River 

0.50 

0.66 

0.91 

1.68 

0.86 

0.65 

0.83 

0.87 

McCormick 

2.85 

2.89 

3.79 

4.18 

1.80 

3.55 

2.69 

3.11 

Oregon 

0.18 

0.15 

0.08 

0.28 

0.05 

0.18 

0.21 

0.16 

Alligator 

0.43 

3.41 

4.74 

4.40 

0.24 

0.62 

0.42 

2.04 

Stillwater 

0.21 

0.08 

0.00 

0.21 

0.00 

0.00 

0.00 

0.07 

East  Highway 

0.11 

0.01 

0.00 

0.00 

0.00 

0.00 

0.00 

0.02 

Embankment 

0.15 

0.15 

0.16 

10.76 

0.00 

0.00 

0.00 

1.60 

NMP  to  P46 

0.04 

0.09 

0.03 

0.12 

0.05 

0.05 

0.11 

0.07 

P46  to  P67 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

Total 

193.97 

241.13 

206.74 

279.45 

190.00 

217.57 

208.27 

219.59 

Outflows 

Name 

Year 

1996  1997  1998  1999  2000  2001  2002 

Average 

Annual 

Evapotranspiration 

-82.15 

-85.05 

-92.09 

-80.45 

-88.84 

-67.41 

-87.43 

-83.35 

Infiltration 

-39.67 

-38.93 

-30.99 

-39.32 

-28.96 

-50.72 

-34.32 

-37.56 

Leakage 

-17.24 

-17.11 

-12.40 

-24.93 

-15.14 

-18.66 

-16.46 

-17.42 

C-lll 

-0.46 

-1.23 

-0.38 

-1.52 

-0.57 

-0.58 

-0.44 

-0.74 

Trout 

-27.23 

-37.44 

-29.72 

-36.74 

-28.68 

-27.50 

-31.81 

-31.30 

Mud 

-6.11 

-7.29 

-8.28 

-11.02 

-7.01 

-10.65 

-8.47 

-8.40 

West  Highway 

-6.01 

-6.04 

-7.94 

-6.22 

-6.07 

-5.25 

-4.98 

-6.07 

East  Creek 

-5.65 

-7.19 

-5.75 

-14.26 

-4.58 

-4.08 

-5.22 

-6.68 

Taylor 

-4.07 

-5.32 

-3.88 

-11.38 

-2.91 

-4.05 

-3.45 

-5.01 

McCormick 

-6.37 

-8.43 

-5.39 

-4.31 

-4.50 

-7.97 

-8.98 

-6.56 

Oregon 

-3.30 

-3.12 

-4.70 

-3.74 

-3.12 

-2.83 

-2.62 

-3.35 

Alligator 

-3.94 

-6.66 

-6.84 

-5.88 

-0.24 

-0.75 

-0.48 

-3.54 

Stillwater 

-2.75 

-3.84 

0.00 

-2.73 

0.00 

0.00 

0.00 

-1.33 

East  Highway 

-1.96 

-1.34 

0.00 

0.00 

0.00 

0.00 

0.00 

-0.47 

Embankment 

-1.40 

-0.05 

-3.50 

-11.83 

0.00 

-0.02 

0.00 

-2.40 

NMP  to  P46 

-0.03 

-0.02 

-0.02 

-0.02 

-0.03 

-0.01 

-0.02 

-0.02 

P46  to  P67 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

Total 

-208.34 

-229.05 

-211.89 

-254.34 

-190.64 

-200.49 

-204.70 

-214.21 
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4.3.1  Surface-Water  Flow  Patterns 

The  average  annual  water  budget  for  the  coastal  wetland  shows  the  relative  magnitude  of  the  different  hydrologic  processes 
and  the  spatial  distribution  of  surface-water  inflows  and  outflows  (fig.  7).  Of  the  average  annual  rainfall  total  (138  cm),  83  cm  is 
lost  to  surface-water  evaporation  and  38  cm  is  lost  to  direct  infiltration  to  the  water  table  as  a  result  of  dry  surface-water  cells.  The 
remaining  17  cm  combines  with  boundary  inflows  from  Taylor  Slough  Bridge  (13.7  cm),  L-31W  (4.3  cm),  and  C-lll  (25.0  cm) 
and  discharges  into  Florida  Bay  through  the  coastal  creeks. 

Simulation  results  clearly  indicate  the  presence  of  spatial  and  temporal  variations  in  wetland  flow  velocity.  Vector  plots  of 
daily  average  flow  velocities  were  constructed  for  November  20,  1999,  and  June  24,  2000,  to  show  wetland  flow  patterns  under 
wet  and  dry  conditions,  respectively  (figs.  8  and  9).  Vector  colors  were  scaled  from  0  to  5  mm/s  to  highlight  wetland  flow  velocity 
only,  rather  than  wetland  and  estuary  flow  velocities,  which  encompass  a  much  larger  range.  Arrows  are  shown  for  all  active  sur¬ 
face-water  cells  with  standing  water.  Figure  8  shows  flow  velocities  for  November  20,  1999 — a  day  with  relatively  high  water 
levels  about  33  days  after  Hurricane  Irene,  which  deposited  more  than  30  cm  of  rain  in  about  46  hours  (Knight  and  others,  2000). 
Taylor  Slough  clearly  is  a  predominant  hydrologic  feature,  with  simulated  flow  velocities  that  exceed  5  mm/s.  Flow  within  Taylor 
Slough  does  not  appear  to  continue  toward  the  southwest  to  the  Cuthbert  and  Seven  Palms  Lake  area,  but  instead  appears  to  be 
routed  toward  the  southeast  into  Taylor  River,  East  Creek  and  Joe  Bay  (fig.  4).  Simulated  flow  velocities  in  the  eastern  part  of  the 
wetland  model  domain  also  exceed  5  mm/s.  These  relatively  large  velocities  are  the  result  of  substantial  inflows  from  the  water- 
management  canals  (L-31W  and  C-lll)  that  were  actively  draining  urban  and  agricultural  areas  in  the  northeast. 

For  June  24,  2000,  a  day  with  relatively  dry  conditions,  standing  surface  water  is  absent  over  much  of  the  area.  Over  the  north¬ 
ern  half  of  the  model  domain,  large  areas  adjacent  to  the  slough  are  dry.  The  topographically  high  areas  surrounding  the  south¬ 
western  lakes  are  also  dry  (fig.  9).  Although  Taylor  Slough  has  standing  water,  flow  velocities  are  less  than  1  mm/s.  On  the  eastern 
side  of  the  model,  there  is  inflow  from  the  C-lll  canal,  but  in  general,  wetland  flow  velocities  for  June  24,  2000,  are  5  to  10  times 
less  than  for  November  20,  1999.  Model  results  also  show  that  simulated  flow  velocities  within  Florida  Bay  generally  are  about 
1  to  3  orders  of  magnitude  greater  than  wetland  flow  velocities  because  of  the  lower  frictional  resistance  values  (Manning’s  n)  and 
larger  depths  assigned  to  Florida  Bay. 


Figure  8.  Daily  average  surface-water  flow  velocities  for  November  20, 1999,  during  a  relatively  wet  period. 
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Figure  9.  Daily  average  surface-water  flow  velocities  for  June  24,  2000,  during  a  relatively  dry  period. 


4.3.2  Discharge  to  Florida  Bay  and  Surface-Water/Ground-Water  Interactions 


Discharge  of  fresh  or  brackish  water  from  the  Florida 
mainland  into  northeastern  Florida  Bay  can  occur  in  three 
ways:  (1)  discharge  from  coastal  creeks,  (2)  overtopping  of 
the  Buttonwood  Embankment,  and  (3)  submarine  ground- 
water  discharge.  Of  these  three  discharge  mechanisms, 
coastal  creek  discharge  is  the  only  one  that  has  been  directly 
measured  in  the  field,  and  continuous  measurements  for 
1996  to  2002  are  available  at  5  of  the  10  coastal  creeks  in  the 
area.  For  the  five  creeks  with  continuous  discharge  measure¬ 
ments  for  the  7-year  simulation  period  (McCormick,  Mud, 
Trout,  Taylor  River,  and  West  Highway),  the  simulated 
cumulative  discharge  is  only  about  15  percent  less  than  the 
measured  cumulative  discharge  (fig.  10).  Based  on  model 
results,  cumulative  discharge  estimates  for  the  five  coastal 
creeks  without  continuous  record  for  the  simulation  period 
(Alligator,  East,  Stillwater,  Oregon,  and  East  Highway) 
comprise  about  24  percent  of  the  measured  cumulative  dis¬ 
charge  at  the  five  monitored  creeks.  The  frictional  resistance 
parameters  of  the  creeks  without  continuous  data  were  not 
altered  from  their  original  values,  which  were  assigned  based 
on  field  observations  of  creek  widths  and  literature  values  for 
roughness  coefficients. 
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Figure  10.  Measured  and  simulated  cumulative  discharge  for 
the  five  measured  creeks.  Trout  Creek,  and  the  five  creeks 
without  continuous  measurements  for  the  7-year  simulation 
period  (1996-2002). 
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Trout 


Figure  11.  Ranking  of  the  coastal  creeks  by  cumulative  dis¬ 
charge  volume  to  Florida  Bay  for  the  period  from  1996  to  2002. 
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Figure  12.  Measured  and  simulated  monthly  discharge  atTrout 
Creek  from  1996  to  2002. 


LU 


O 

> 

LU 

0 

GC 

< 

X 

o 

0 

o 

> 


< 

0 


J  FMAMJ  J  ASOND 


1999 


Figure  13.  Measured  and  simulated  average  daily  discharge 
atTrout  Creek  for  1999. 


The  distribution  of  coastal  creek  discharge  to  Florida  Bay  was 
evaluated  by  using  measured  discharge  volumes  for  the  five  creeks 
with  continuous  data  for  the  7-year  simulation  period  and  simulated 
discharge  volumes  for  the  remaining  five  creeks.  For  the  7-year 
simulation  period,  Trout  Creek  contributed  nearly  half  (47  percent) 
of  the  discharge  to  Florida  Bay  (fig.  11).  West  Highway,  East,  Tay¬ 
lor  River,  and  Mud  Creeks  were  the  next  largest  contributors  with 
12, 10,  8  and  7  percent  of  the  total  discharge,  respectively.  The  five 
remaining  creeks  (Oregon,  Alligator,  Stillwater,  McCormick,  and 
East  Highway)  each  contributed  5  percent  or  less  of  the  total  dis¬ 
charge. 

The  coastal  creeks  show  distinct  seasonal  patterns  in  discharge 
to  Florida  Bay.  For  example,  at  the  monthly  timescale,  simulated 
and  measured  discharges  at  Trout  Creek  peak  during  the  wet  season 
and  reverse  directions  during  most  dry  seasons  (fig.  12).  Discharge 
volumes  begin  to  increase  with  the  onset  of  the  wet  season,  which 
typically  occurs  between  May  and  June.  By  September  or  October 
of  each  year,  discharge  volumes  reach  their  annual  peaks  and  then 
begin  to  decline  as  the  dry  season  approaches  in  November.  During 
some  dry  season  months,  discharge  rates  are  negative,  which  indi¬ 
cate  northward  flow  from  Florida  Bay  into  the  coastal  wetlands. 
These  negative  discharge  rates,  which  are  reproduced  by  the  model, 
are  primarily  caused  by  sustained  periods  of  southerly  winds  that 
push  brackish  Florida  Bay  water  inland.  The  ability  of  the  model  to 
match  these  negative  discharge  values  proved  to  be  critical  in  accu¬ 
rately  representing  salinity  values  in  the  coastal  wetlands.  For  the 
entire  7-year  period,  the  average  monthly  discharge  volume  calcu¬ 
lated  with  measured  data  is  1.7  x  107  m3  (fig.  12).  The  average 
monthly  discharge  volume  calculated  using  simulated  data  is  1.2  x 
107  m3,  about  27  percent  too  low.  Using  the  measured  and  simu¬ 
lated  records  in  fig.  12,  the  R2  value  is  0.77.  This  discrepancy 
appears  to  be  caused  by  failure  of  the  model  to  capture  peak  flows 
during  the  wet  season.  A  possible  explanation  is  that  one  or  more 
of  the  inflows  represented  by  the  model  is  based  on  inaccurate  field 
data.  The  model  is  highly  sensitive  to  rainfall,  and  small  errors  in 
rainfall  values,  when  applied  to  a  large  area,  can  lead  to  substantial 
errors  in  creek  discharge.  At  shorter  timescales,  negative  discharge 
values  occur  frequently  throughout  the  year.  Figure  13  shows  daily 
discharge  volumes  at  Trout  Creek  for  1999.  Again,  the  model 
seems  to  represent  the  range  in  discharge  volumes,  capturing  both 
the  high  and  low  values.  The  model  does,  however,  fail  to  represent 
some  of  the  higher  peaks,  which  results  in  a  20-percent  underesti¬ 
mation  of  average  annual  discharge  volume  at  Trout  Creek  for 
1999.  Using  daily  average  flows,  the  R2  value  for  measured  and 
simulated  discharges  at  Trout  Creek  is  0.78. 

Spectral  analysis  was  performed  on  discharge  data  at  Trout 
Creek  to  determine  if  the  dominant  frequencies  observed  in  the  field 
measurements  are  represented  by  the  model.  Figure  14  shows  an 
amplitude  spectrum  of  measured  and  simulated  discharge  for  the 
7-year  dataset  at  Trout  Creek.  Four  distinct  spikes  in  the  spectrum 
are  shown  at  frequencies  of0.93, 1.00, 1.93,and2.00  cycles  per  day 
in  both  the  measured  and  simulated  discharges.  The  spike  at  0.93 
corresponds  to  the  01  tidal  component,  which  has  a  period  of 
25.8  hours.  At  this  frequency,  the  amplitude  spectrum  values  for 
the  measured  and  simulated  discharges  are  0.209  and  0.208  m3/d, 
respectively.  The  1.00  and  2.00  cycles  per  day  frequencies 
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correspond  to  periods  of  24  and  12  hours,  respectively,  and  are  caused  by  temporal  variations  in  wind,  and  possibly  the  SI  and 
S2  solar  tides.  Dominant  spikes  at  the  1.00  and  2.00  frequencies  are  also  seen  in  the  spectrum  of  wind  velocity  (not  shown)  col¬ 
lected  at  the  Joe  Bay  weather  station.  At  the  1.00  cycle  per  day  frequency,  the  amplitude  spectrum  values  for  the  measured  and 
simulated  discharge  are  0.529  and  0.763  m3/d,  respectively.  At  the  2.00  cycle  per  day  frequency,  the  amplitude  spectrum  values 
for  the  measured  and  simulated  discharge  are  0.394  and  0.428  m3/d,  respectively.  The  spike  at  the  1.93  frequency  corresponds  to 
the  M2  lunar  tide.  At  this  frequency,  the  simulated  amplitude  spectrum  (0.199  m3/d)  is  less  than  half  of  the  measured  amplitude 
spectrum  (0.5 12  m3/d).  A  harmonic  analysis  of  Trout  Creek  stage  indicates  that  the  M2  amplitude  is  0.38  cm,  which  is  only  slightly 
larger  than  the  precision  of  the  stage  recorder  (0.30  cm).  Thus,  a  possible  explanation  for  the  discrepancy  between  the  simulated 
and  measured  amplitude  spectrum  at  the  M2  frequency  is  that  the  boundary  stages  in  the  model  are  not  recorded  with  enough  pre¬ 
cision  to  reproduce  the  complete  M2  signal. 

Buttonwood  Embankment  overtopping  is 
another  mechanism  that  discharges  freshwater 
from  the  coastal  wetlands  into  Florida  Bay.  Due 
to  the  remote  location  and  expansive  length  of 
the  embankment,  however,  overtopping  dis¬ 
charge  volumes  have  never  been  measured.  Ele¬ 
vations  of  the  embankment  crest  are  not  avail¬ 
able,  except  at  the  coastal  creeks  where  estimates 
can  be  made  based  on  the  height  of  the  embank¬ 
ment  above  the  water  surface.  In  the  model  grid, 
the  embankment  height  is  set  at  0.3  m  above  land 
surface.  Simulation  results  indicate  that  embank¬ 
ment  overtopping  is  infrequent,  but  can  occur  in 
both  directions  in  response  to  tropical  storms. 

For  example,  as  Hurricane  Irene  approached  the 
Florida  mainland  in  October  1999,  a  storm  surge 
was  recorded  at  the  Taylor  River  monitoring  sta¬ 
tion  (fig.  15).  This  storm  surge  pushed  a  large 
volume  of  brackish  water  from  Florida  Bay  over 
the  embankment  and  into  the  coastal  wetlands. 

After  Hurricane  Irene  made  landfall,  the  associ¬ 
ated  heavy  rainfall  reversed  flow  over  the 
embankment  and  into  Florida  Bay,  as  indicated 
by  the  positive  discharge  values  (fig.  15).  This 
was  the  largest  overtopping  event  simulated  by 
the  model.  For  the  entire  7 -year  simulation 
period,  the  net  overtopping  discharge  volume 
was  -3.7  x  107  m3.  This  suggests  that  although 
overtopping  may  allow  for  flow  into  the  coastal 
wetlands,  the  mechanism  is  not  a  substantial 
source  of  freshwater  for  Florida  Bay.  The  cumu¬ 
lative  positive  and  negative  overtopping  dis¬ 
charge  volumes  are  7.3  x  107  m3  and  -1.1  x  108 
m3,  respectively,  which  are  relatively  small  com¬ 
pared  to  the  cumulative  creek  volumes  for  the  7- 
year  simulation  period  (about  28  x  108  m3  for  the 
10  coastal  creeks).  At  daily,  weekly,  or  monthly 
timescales,  however,  the  overtopping  volumes 
may  be  significant  in  terms  of  freshwater  flows 
into  Florida  Bay,  or  brackish  water  flow  into  the 
coastal  wetland. 

Daily  leakage  rates  between  surface  water  and  ground  water  are  produced  as  part  of  the  model  output  for  each  cell.  These 
daily  leakage  rates  were  averaged  over  the  7-year  simulation  period  to  illustrate  the  spatial  variability  of  surface-water/ground- 
water  interaction  and  to  determine  whether  ground  water  is  discharging  into  Florida  Bay.  These  leakage  rates  do  not  include 
recharge  or  evapotranspiration  directly  to  or  from  the  water  table.  The  model  suggests  an  alternating  pattern  of  downward  and 
upward  leakage  from  north  to  south  (fig.  16).  Within  the  wetland  portion  of  the  model  domain,  downward  leakage  rates  (shown  as 


Figure  14.  Amplitude  spectrum  as  a  function  of  frequency  calculated  using  2-hour 
measured  and  simulated  discharges  at  Trout  Creek  for  the  7-year  simulation  peri¬ 
od. 
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Figure  15.  Discharge  over  Buttonwood  Embankment  and  stage  at  Taylor 
River  during  Hurricane  Irene,  October  1999. 
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Figure  16.  Average  annual  leakage  rates  for  the  7-year  simulation  period. 


positive  values)  correlate  with  topographically  high  areas.  For  example,  in  the  northern  part  of  the  model,  average  leakage  rates 
exceed  0.2  cm/d.  A  band  of  upward  leakage  with  rates  in  some  areas  exceeding  0.2  cm/d  appears  across  the  central  part  of  the 
model.  This  upward  leakage  band  correlates  with  the  topographically  low  area  within  the  central  part  of  the  model  (fig.  5)  and  with 
the  freshwater/saltwater  interface  in  the  aquifer  (fig.  4).  Downward  leakage  rates  also  occur  just  north  of  the  Buttonwood  Embank¬ 
ment  where  surface-water  levels  tend  to  be  higher  than  those  in  Florida  Bay.  These  model  results  suggest  that  there  may  be  shallow 
ground- water  flow  beneath  the  embankment,  which  then  discharges  into  Florida  Bay.  The  source  for  this  shallow  ground- water 
flow  system  is  surface  water  impounded  by  the  Buttonwood  Embankment.  South  of  the  Buttonwood  Embankment,  ground  water 
discharges  upward  into  the  coastal  embayments  of  Florida  Bay.  Average  leakage  rates  within  this  zone  can  exceed  0.2  cm/d,  but 
most  are  between  0.01  and  0.1  cm/d.  Model  results  also  indicate  that  Joe  Bay  (figs.  4  and  16)  may  be  a  ground- water  discharge 
area.  Joe  Bay  has  an  exposed  limestone  bottom,  and  thus,  the  absence  of  a  thin  layer  of  low-permeability  sediments  results  in  a 
relatively  strong  hydraulic  connection  between  surface  water  and  the  underlying  aquifer.  Average  leakage  rates  appear  to  be  down¬ 
ward  over  most  of  Florida  Bay  (figs.  4  and  16),  with  values  ranging  between  about  0.0  and  0.2  cm/d.  Downward  leakage  in  this 
zone  is  probably  the  result  of  cyclic  flow  that  often  occurs  in  freshwater/saltwater  interfaces  within  a  coastal  aquifer  (Kohout,  1964; 
Langevin,  2001).  Fresh  ground  water  flowing  toward  an  interface  mixes  with  saline  ground  water.  This  brackish  mixture  then 
discharges  into  the  ocean,  coastal  estuary,  or  in  this  case,  into  the  brackish  water  wetlands. 

Average  leakage  rates  for  the  entire  simulation  period  are  shown  in  (fig.  16);  however,  daily  leakage  rates  are  highly  variable 
and  can  change  direction  in  response  to  rainfall  events  or  prolonged  dry  periods.  Mapped  leakage  patterns  for  specific  days  and 
averaged  periods  are  similar,  except  after  large  rainfall  events.  Graphs  of  simulated  daily  leakage  and  water  levels  were  prepared 
for  selected  model  cells  (labeled  A-D  in  fig.  16).  The  average  leakage  rate  is  about  -0.05  cm/d  at  point  A  (fig.  17),  which  is  in  a 
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Figure  17.  Simulated  leakage  rates,  surface-water  stage,  and  ground-water  head  at  four  locations. 
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zone  of  upward  leakage  just  south  of  Old  Ingraham  Highway  (fig.  16).  During  most  of  the  simulation  period,  leakage  at  point 
A  is  upward.  Large  rainfall  events,  such  as  in  June  1997,  however,  seem  to  greatly  affect  the  vertical  movement  of  water.  Ver¬ 
tical  leakage  in  the  area  during  these  events  appears  to  change  direction,  with  surface-water  flowing  downward  into  the  aquifer. 
Temporally  variable  leakage  patterns  similar  to  location  A  also  are  shown  for  points  B  and  C  (fig.  17);  downward  leakage  rates 
into  the  aquifer  are  relatively  high  after  large  rainfall  events.  During  periods  when  surface-water  levels  decline  between  rainfall 
events,  ground  water  appears  to  leak  gradually  upward  and  mix  with  the  surface  water.  Point  D  in  Florida  Bay  (fig.  16)  appears 
to  respond  differently  than  the  other  locations,  possibly  because  it  is  highly  affected  by  the  stage  and  salinity  of  Florida  Bay 
(fig.  17).  Simulated  ground- water  leakage  rates  at  location  D  range  from  -1.7  to  1.7  cm/d,  with  an  average  rate  of  -0.09  cm/d. 
The  water  budget  for  the  coastal  wetland  part  of  the  model  domain  indicates  that  average  volumes  of  ground-water  recharge 
and  discharge  are  nearly  identical  for  the  simulation  period  (fig.  7  and  table  3).  For  any  particular  year,  however,  the  wetland 
may  experience  a  net  loss  or  gain  as  a  result  of  leakage  (table  3). 

Simulated  exchange  rates  between  surface  water  and  ground  water  contain  some  degree  of  uncertainty,  with  perhaps  the 
exception  of  northern  Taylor  Slough,  because  the  model  was  not  calibrated  to  direct  leakage  measurements.  Harvey  and  others 
(2000a)  suggest  that  leakage  rates  into  the  northern  part  of  Taylor  Slough  (south  of  Old  Ingraham  Highway)  may  be  as  much 
as  3  cm/d.  The  simulated  leakage  rate  (about  0.25  cm/d  for  the  same  area  and  month),  however,  is  an  order  of  magnitude  lower 
than  the  maximum  measured  rate,  but  is  in  the  correction  direction  (upward)  and  probably  within  the  error  range  of  the  chloride 
dilution  method  used  for  the  calculation.  Harvey  and  others  (2000b)  also  measured  head  differences  across  the  peat  layer  at 
1 1  locations  during  six  different  field  visits.  To  compare  these  field  measurements  with  the  model,  simulated  ground- water 
heads  from  layer  1  were  converted  from  equivalent  freshwater  heads  to  aquifer  heads  and  then  subtracted  from  overlying  sur¬ 
face-water  stages  to  calculate  simulated  head  differences  across  the  peat.  Numerous  simulations  were  performed  with  varying 
hydraulic  properties  for  the  peat  and  uppermost  aquifer  layer  in  an  effort  to  calibrate  the  model.  For  the  36  head-difference 
measurements,  the  mean  and  mean  absolute  errors  between  the  simulated  and  measured  values  are  -0.9  and  5.2  cm,  respec¬ 
tively.  Comparisons  of  simulated  head  differences  with  the  observed  head  differences  indicate  that  many  of  the  simulated 
directions  of  vertical  leakage  are  correct.  Discrepancies  between  simulated  and  measured  values  are  likely  caused  by  variations 
in  field  head  differences  that  occur  over  distances  shorter  than  the  305-m  cell  size  used  in  the  model,  or  by  the  presence  of  a 
complex  and  heterogeneous  peat  layer  with  spatially  variable  hydraulic  properties. 

4.3.3  Coastal  Salinities 

Field  data  and  model  results  indicate  a  strong  seasonal  pattern  in  coastal  wetland  salinities.  Figures  18  and  19  are  maps  of 
daily  average  surface-water  salinities  for  November  20,  1999,  and  for  June  24,  2000,  respectively.  These  dates  are  the  same 
as  used  in  figures  8  and  9  to  show  daily  average  surface-water  flow  velocities.  During  the  wet  season,  salinities  for  most  of  the 
coastal  wetlands  are  nearly  zero,  except  for  several  isolated  areas  north  of  the  Buttonwood  Embankment  where  salinities  are 
less  than  about  10  psu  (figs.  4  and  18).  In  Florida  Bay,  salinities  on  November  20,  1999,  are  much  less  than  the  salinity  of 
seawater.  For  example,  in  Little  Madeira  Bay,  salinities  near  the  mouths  of  East  Creek  and  Taylor  River  are  about  5  to  10  psu. 
Salinities  farther  offshore  range  between  about  15  and  25  psu.  The  salinity  pattern  is  quite  different  for  June  24,  2000  (fig.  19). 
Salinities  from  5  to  10  psu  have  extended  into  the  coastal  wetlands  west  of  Joe  Bay.  The  simulated  salinity  in  Joe  Bay  ranges 
between  about  15  and  23  psu  (figs.  4  and  19).  Most  of  the  Florida  Bay  area,  except  for  Little  Madeira  Bay,  Long  Sound,  and 
Little  Blackwater  Sound  (shown  on  fig.  4),  has  seawater  salinities  of  about  35  psu.  A  comparison  between  measured  and  sim¬ 
ulated  values  of  average  monthly  salinity  at  Trout  Creek  is  shown  in  figure  20.  The  field  measurements  and  simulation  clearly 
indicate  season  fluctuations  in  salinity  at  Trout  Creek.  Salinities  reach  35  psu  in  May,  June,  or  July  of  each  year,  which  corre¬ 
sponds  with  the  end  of  the  dry  season.  The  lowest  salinities  were  recorded  in  August,  September,  or  October  of  each  year. 
Using  the  average  monthly  data  in  figure  20,  the  measured  and  simulated  salinities  have  an  R2  value  of  0.76.  A  comparison  of 
daily  salinities  gives  an  R2  value  of  0.67,  suggesting  that  the  model  is  better  at  representing  the  longer  seasonal  fluctuations 
than  the  shorter  timescale  daily  or  weekly  salinity  fluctuations. 
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Figure  18.  Daily  average  surface-water  salinities  for  November  20, 1999,  during  a  relatively  wet  period. 


Figure  19.  Daily  average  surface-water  salinities  for  Junen  24,  2000,  during  a  relatively  dry  period. 
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Figure  20.  Measured  and  simulated  values  of  monthly  average 
salinity  at  Trout  Creek  for  the  7-year  simulation  period  (1996-2002) 


4.4  Effects  of  Selected  Hydrologic  Processes 


Several  simulations  were  performed  to  evaluate  the  effects  of  hydrologic 
processes  unique  to  this  particular  model  application,  namely  (1)  surface-water 
and  ground- water  interactions,  (2)  density-dependent  flow,  and  (3)  local  wind 
stress.  These  three  processes  are  all  known  to  be  active  within  the  study  area, 
and  thus,  their  effects  can  be  evaluated  by  comparing  simulations  without  these 
processes  to  the  previously  described  integrated  simulation,  referred  to  here  as 
the  base  case.  A  possible  limitation  with  this  approach  is  that  the  consequences 
for  neglecting  a  process  may  be  overstated  if  calibration  of  the  integrated  model 
tended  to  overemphasize  that  process. 


Table  4.  Calibration  statistics  for  the  surface- 
water  simulation  without  leakage.  Errors  are  cal¬ 
culated  relative  to  field  data  for  daily  average 
coastal  creek  discharges  and  coastal  creek  sa¬ 
linities.  Mean  errors  were  calculated  by  sub¬ 
tracting  measured  values  from  simulated  values. 
Station  locations  are  shown  in  figure  4. 

[ME,  mean  error;  MAE,  mean  absolute  error;  RMSE, 
root  mean  square  error;  m3/s,  cubic  meters  per  second;  m, 
meter;  psu,  practical  salinity  units] . 


The  ground- water  part  of  the  integrated  model  was  calibrated  to  heads  at 
three  monitoring  wells  (G-3619,  G-3353,  G-1251),  36  head-difference  measure¬ 
ments,  the  position  of  the  freshwater/saltwater  interface,  and  to  estimates  of 
leakage.  Further  adjustments  to  the  ground- water  model  did  not  improve  simu¬ 
lation  results,  and  thus  the  integrated  model  was  considered  calibrated  within  the 
limitations  of  the  trial  and  error  method.  The  effect  of  leakage  on  the  surface- 
water  system  was  evaluated  by  removing  the  ground- water  model.  Neglecting 
surface-water  and  ground-water  interactions  tends  to  worsen  simulated  dis¬ 
charges  and  salinities  in  most  cases.  Mean  absolute  errors  (MAE)  and  root  mean 
square  errors  (RMSE)  for  the  simulation  without  leakage  (table  4)  are  larger 
(with  the  exception  of  salinity  errors  at  McCormick  Creek  and  Taylor  River) 
than  errors  for  the  base  case  (table  2).  The  average  increases  in  MAE  and  RMSE 
for  coastal  creek  discharge  as  a  result  of  neglecting  leakage  are  0.54  and 

Q 

0.94  m  /d,  respectively.  The  average  increases  in  MAE  and  RMSE  for  coastal 
creek  salinity  are  0.14  and  0.19  psu,  respectively. 

In  the  second  simulation,  fluid  density  was  held  constant  in  space  and  time  by  adjusting  the  equation  of  state  (eq.  5)  such  that 
concentration  did  not  affect  fluid  density.  The  resulting  cumulative  flow  through  the  five  measured  creeks  (1.73  x  109  m3)  is  about 
9  percent  less  than  for  the  base  case  simulation  (1.91  x  109  m3)  and  about  24  percent  less  than  the  measured  cumulative  discharge 
(2.25  x  10  m  ).  There  is  little  difference  in  cumulative  discharge  for  the  five  creeks  without  continuous  data  (4.60  x  10  m  com¬ 
pared  with  5.43  x  108  m3  for  the  base  case).  Some  minor  differences  between  this  simulation  and  the  base  case  were  noted  for 
leakage  rates,  but  in  general  the  leakage  pattern  is  similar  to  that  for  the  base  case  (fig.  16).  These  results,  therefore,  indicate  that 
the  upward  leakage  zone  located  between  A  and  C  in  figure  16  is  caused  by  topographic  variations  rather  than  by  variable-density 
effects  near  the  relatively  dense  saltwater  wedge  observed  in  the  Biscayne  aquifer. 


Station 

ME 

MAE 

RMSE 

Discharge  (m3/s) 

McCormick 

0.78 

2.62 

4.98 

Mud 

0.25 

2.40 

4.32 

Trout 

-2.13 

5.74 

7.86 

Taylor  River 

-0.06 

1.50 

3.40 

West  Highway 

-0.29 

1.18 

1.65 

Salinity  (psu) 

McCormick 

2.19 

6.18 

7.53 

Mud 

-2.05 

5.04 

6.48 

Trout 

2.36 

5.45 

7.00 

Taylor  River 

1.78 

4.84 

5.74 

West  Highway 

4.97 

6.10 

8.45 
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The  local  wind  stress  was  removed  for  the  third  simulation. 
A  distinction  is  made  here  between  local  and  regional  wind 
effects.  In  SWIFT2D,  the  local  wind  stresses  are  included  in  the 
conservation  of  momentum  equations  (eqs.  2  and  3).  The  model 
allows  input  of  temporally  and/or  spatially  varying  wind  speed 
and  direction  for  calculation  of  stress.  On  the  other  hand, 
regional  wind  effects  are  included  in  the  limited  domain  model 
through  specified  water-level  boundaries.  For  example,  a  strong 
southerly  wind  over  Florida  Bay  will  push  water  against  the  But¬ 
tonwood  Embankment  and  raise  water  levels  in  northeastern 
Florida  Bay.  Thus,  the  water  levels  measured  in  northeastern 
Florida  Bay,  which  are  used  to  assign  the  southernmost  stage 
boundary  for  the  model,  contain  the  effect  of  the  regional  wind 
field.  Removing  the  local  wind  stress  does  not  have  a  substantial 
effect  on  coastal  creek  flows,  but  does  affect  coastal  salinities. 
As  an  example,  daily  average  salinities  at  Trout  Creek  from 
November  2001  to  July  2002  are  shown  in  figure  21.  Clearly, 
the  simulation  is  improved  when  the  local  wind  stress  is 
included  in  the  model,  both  in  terms  of  the  short-term  fluctua¬ 
tions  observed  at  the  end  of  the  2001  wet  season  and  in  terms  of 
the  longer  time  increase  in  salinity  as  a  result  of  the  dry  season. 


DATE 

Figure  21.  Average  daily  salinity  atTrout  Creekfrom  simulation 
without  local  wind  stress,  the  base  case  simulation,  and  from 
measured  data. 


5.  Discussion 

Prior  to  performing  simulations  with  the  integrated  model,  the  surface-water  and  ground- water  models  were  independently 
developed  and  calibrated  to  the  extent  possible.  For  the  initial  surface-water  model,  exchange  with  ground  water  was  considered 
negligible  (Swain  and  others,  2004).  Ground- water  model  development  was  performed  after  the  surface-water  model  was  devel¬ 
oped,  and  thus  simulated  surface-water  stages  and  salinities  were  applied  as  boundary  conditions  over  the  aquifer  surface.  This 
stepwise  approach  had  two  advantages.  First,  it  was  relatively  easy  to  identify  and  correct  input  and  runtime  errors  for  the  individ¬ 
ual  models  before  they  were  integrated.  Second,  computer  runtimes  for  the  ground- water  model  were  only  a  couple  of  hours, 
whereas  the  integrated  model  required  over  30  hours  to  run.  The  shorter  computer  runtimes  were  particularly  useful  during  cali¬ 
bration  of  the  ground-water  model  to  aquifer  salinity.  Because  of  the  highly  transmissive  nature  of  the  Biscayne  aquifer  and  a  rel¬ 
atively  stable  freshwater/saltwater  interface  in  southern  Florida  (Sonenshein,  1997),  aquifer  salinities  were  assumed  to  be  in  equi¬ 
librium  with  current  water  levels  and  hydrologic  stresses.  Thus,  an  additional  level  of  confidence  in  the  ground- water  model  was 
established  when  it  could  be  shown  that  after  the  model  reached  dynamic  equilibrium  (through  repeated  simulations),  the  simulated 
freshwater/saltwater  interface  was  in  the  observed  location.  Only  minor  salinity  adjustments  at  the  ground- water  boundaries  were 
required  as  part  of  this  calibration  process  as  hydrodynamic  dispersion  was  not  active  for  the  simulations. 

Limitations  were  periodically  encountered  using  the  explicit,  time-lagged  approach  to  couple  the  surface-water  and  ground- 
water  models.  For  some  sensitivity  simulations  with  very  large  leakage  rates,  convergence  could  not  be  achieved  during  solution 
of  the  ground- water  flow  equation.  Evaluation  of  the  convergence  problems  indicated  that  very  large  leakage  rates  caused  numer¬ 
ical  oscillations  in  the  implicit  solution.  Ground-water  heads  measured  in  the  field  can  respond  quickly  to  hydrologic  stresses.  For 
the  Everglades  application  of  the  integrated  model,  however,  large  leakage  rates  may  persist  throughout  the  day  in  the  model, 
because  of  the  1-day  length  of  the  stress  period  in  SEAWAT,  whereas  actual  leakage  rates  would  decrease  as  ground-water  heads 
respond  more  quickly.  These  convergence  problems  could  probably  have  been  avoided  by  decreasing  the  length  of  the  ground- 
water  stress  period.  These  convergence  problems  were  encountered  only  in  a  few  instances,  and  thus  the  day  lag,  which  is  compu¬ 
tationally  many  times  faster  than  using  an  hourly  lag  or  fully  implicit  solution,  was  a  necessity  for  this  particular  application.  Future 
efforts  using  the  integrated  model  should  follow  the  example  of  Fairbanks  and  others  (2001)  and  focus  on  determining  the  relation 
between  accuracy  and  efficiency  for  different  coupling  approaches  and  timestep  lengths. 

The  Buttonwood  Embankment  clearly  is  an  important  physiographic  feature  in  the  Taylor  Slough  area.  Model  results  and  field 
observations  suggest  that  freshwater  flow  into  Florida  Bay  occurs  primarily  through  the  coastal  creeks,  rather  than  as  overtopping 
of  the  embankment.  This  flow  pattern  has  allowed  field  investigations  to  quantify  with  a  high  level  of  certainty  the  flow  exchanges 
between  the  coastal  wetlands  and  Florida  Bay  (Hittle,  2000;  Hittle  and  others,  2001).  Confidence  in  the  predictive  capability  of  the 
integrated  model  is  due  largely  to  the  accuracy  and  long-term  record  length  of  creek  discharge  data.  Many  coastal  wetlands  in  other 
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locations,  however,  are  not  separated  from  the  adjacent  marine  water  body  by  an  embankment  or  barrier  to  overland  flow,  and  thus 
measurement  of  freshwater  outflows  is  not  as  straightforward. 

The  integrated  model  described  here  represents  numerous  hydrologic  processes  and  requires  extensive  data  input  to  simulate 
the  flow  and  salinity  patterns  for  the  7-year  period.  Consequently,  the  model  is  subject  to  numerous  limitations,  and  model  results 
should  be  used  with  caution,  particularly  results  for  which  there  are  no  measured  values.  For  example,  estimates  of  Buttonwood 
Embankment  overtopping  discharge  volumes  contain  a  high  level  of  uncertainty,  because  the  crest  elevation  of  the  embankment  is 
more  variable  than  specified  in  the  model.  Simulated  leakage  rates  also  contain  a  high  level  of  uncertainty,  particularly  those  in 
the  southern  part  of  the  coastal  wetland  and  northeastern  Florida  Bay  where  field  estimates  of  surface-water  and  ground- water  inter¬ 
actions  are  lacking.  Future  applications  that  use  forcing  conditions  outside  the  range  of  values  used  for  calibration  would  also  be 
subject  to  limitations.  Although  these  limitations  are  present,  comparisons  between  simulated  and  observed  flow  and  salinity  pat¬ 
terns  in  both  the  wetland  and  aquifer  indicate  that  important  system  processes  and  behavior  are  represented  by  the  model.  The 
model  clearly  reproduces  the  gradual  salinity  increase  during  the  dry  season  and  subsequent  freshwater  flush  during  the  wet  season. 
For  the  7-year  simulation  period,  simulated  rates  of  freshwater  discharge  to  Florida  Bay  are  only  15  percent  less  than  measured 
discharges.  For  these  reasons,  the  model  seems  well  suited  to  serve  its  intended  purpose  of  predicting  the  effects  of  Everglades 
restoration  on  flows,  stages,  and  salinities  in  the  coastal  wetlands,  provided  the  model’s  limitations  are  carefully  considered  in  the 
evaluation  of  predictions. 


6.  Summary  and  Conclusions 

This  paper  describes  a  numerical  approach  for  simulating  integrated  surface-water/ground- water  flow  and  solute  transport  in 
coastal  wetlands  and  adjacent  estuaries.  The  approach  combines  the  SWIFT2D  two-dimensional  hydrodynamic  flow  and  solute- 
transport  code  with  the  SEAWAT  three-dimensional,  saturated  ground-water  flow  and  solute-transport  code.  The  surface-water 
and  ground- water  models,  which  both  simulate  density-dependent  flow,  are  coupled  using  an  explicit  time-lagged  approach  based 
on  a  variable-density  form  of  Darcy’s  law  to  calculate  the  leakage  flux  at  the  ground  surface;  solute  mass  transfer  between  surface 
water  and  ground  water  is  assumed  to  occur  only  by  leakage  advection. 

The  integrated  code  was  applied  to  the  southern  Everglades  of  Florida  and  northeastern  Florida  Bay  to  quantify  flow  and  salin¬ 
ity  patterns  for  the  period  1996-2002  and  to  evaluate  the  effects  of  selected  hydrologic  processes.  The  model  was  calibrated  to  a 
wide  range  of  field  data,  including  coastal  creek  flows  through  the  Buttonwood  Embankment,  a  narrow  but  continuous  feature  that 
separates  Florida  Bay  from  the  coastal  wetlands.  Simulated  surface-water  flow  patterns  in  Taylor  Slough  indicate  southwesterly 
flow  in  the  northern  part  of  the  model  area  with  a  gradual  change  in  flow  direction  toward  the  southeast  into  Joe  Bay  and  through 
Trout  Creek.  For  the  five  creeks  with  continuous  discharge  measurements  for  the  7-year  simulation  period  (McCormick,  Mud, 
Trout,  Taylor  River,  and  West  Highway),  the  simulated  cumulative  discharge  is  only  about  15  percent  less  than  the  measured  cumu¬ 
lative  discharge.  Of  the  10  coastal  creeks,  Trout  Creek  is  the  largest  contributor  of  freshwater  to  Florida  Bay  (47  percent)  with 
West  Highway  (12  percent)  and  East  Creek  (10  percent)  as  the  next  largest  contributors. 

In  addition  to  creek  flows,  the  model  also  simulates  overtopping  of  the  Buttonwood  Embankment  and  submarine  ground- water 
discharge  as  mechanisms  for  delivering  freshwater  from  the  coastal  wetlands  into  Florida  Bay.  Although  simulated  estimates  of 
Buttonwood  Embankment  overtopping  contain  a  high  level  of  uncertainty,  model  results  indicate  that  overtopping  is  infrequent, 
but  can  occur  in  response  to  tropical  storms.  Storm  surges  force  brackish  Florida  Bay  water  over  the  embankment  and  into  the 
coastal  wetlands.  After  making  landfall,  a  tropical  storm  can  also  produce  enough  rain  to  reverse  embankment  overflow  from  the 
coastal  wetland  into  Florida  Bay.  For  the  7-year  simulation  period,  the  net  embankment  overflow,  which  is  into  Florida  Bay,  is 
only  about  1.5  percent  of  the  combined  coastal  creek  flow.  The  water  budget  for  the  coastal  wetland  part  of  the  model  domain 
indicates  that  average  rates  of  downward  leakage  (17.42  cm/yr)  and  upward  leakage  (17.14  cm/yr)  are  nearly  identical  for  the  sim¬ 
ulation  period,  but  for  any  particular  year,  however,  the  wetland  may  experience  a  net  loss  or  gain  to  or  from  the  aquifer.  Model 
results  also  indicate  that  submarine  ground-water  discharge  may  be  occurring  on  the  south  side  of  the  embankment  in  response  to 
the  higher  surface-water  levels  in  the  coastal  wetland.  A  field  survey  would  be  necessary  to  determine  the  validity  of  this  model 
result. 

Field  data  and  model  results  indicate  a  strong  seasonal  pattern  in  coastal  wetland  salinities.  Salinities  at  the  coastal  creeks 
reach  35  psu  toward  the  end  of  the  dry  season,  but  quickly  drop  to  less  than  5  psu  with  the  onset  of  the  wet  season.  This  seasonal 
flushing  pattern  is  well  represented  by  the  model  with  mean  absolute  errors  in  simulated  salinity  ranging  between  4  and  7  psu  for 
the  five  coastal  creeks  with  continuous  data  for  the  7-year  simulation  period.  Future  modifications  to  the  water-management  system 
in  southern  Florida  may  alter  the  freshwater  deliveries  to  the  Taylor  Slough  area.  Based  on  the  performance  of  the  model  to  match 
the  seasonal  flushing  pattern,  the  model  should  be  able  to  predict  the  effects  of  these  altered  water  deliveries  on  coastal  salinity 
patterns. 
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The  effects  of  surface-water  and  ground- water  interactions,  density-dependent  flow,  and  local  wind  stress  were  evaluated  by 
performing  simulations  without  these  processes  and  comparing  results  with  the  base  case  simulation.  In  general,  the  surface-water 
model  that  neglects  interactions  with  ground  water  compares  worse  with  field  data  than  the  base  case  integrated  model;  however, 
without  additional  leakage  measurements,  the  better  match  with  the  integrated  model  cannot  be  conclusively  attributed  to  ground- 
water  interactions.  A  constant-density  simulation  results  in  cumulative  creek  flows  that  are  about  9  percent  less  than  the  base  case, 
and  only  a  slightly  different  pattern  in  leakage,  suggesting  that  the  upward  leakage  zone  that  coincides  with  the  freshwater/saltwater 
interface  in  the  Biscayne  aquifer  is  caused  by  topographic  variations  rather  than  by  density  variations.  Removing  the  local  wind 
stress  does  not  have  a  substantial  effect  on  creek  flows,  but  does  affect  coastal  salinities.  Without  the  local  wind  stress,  Trout  Creek 
salinities  do  not  increase  to  the  30-35  psu  values  measured  in  the  field  during  the  dry  season. 

In  general,  comparisons  between  simulated  and  observed  flow  and  salinity  patterns  in  both  the  wetland  and  aquifer  indicate 
that  important  system  processes  and  behavior  are  represented  by  the  model,  and  although  the  model  is  subject  to  limitations,  it  is 
well  suited  to  predict  the  effects  of  Everglades  restoration  on  the  Taylor  Slough  coastal  wetlands.  The  general  approach  described 
here  would  also  be  applicable  to  other  coastal  wetlands  where  restoration  or  contaminant  transport  issues  are  of  concern.  The  inte¬ 
grated  code  is  robust,  accurate,  and  can  represent  hydrodynamic  surface-water  flow  and  variable-density  ground- water  flow  for 
multi-year  periods. 
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